Implicit type constitutive relations for elastic solids and their use in the development of mathematical models for viscoelastic fluids

Viscoelastic fluids are non-Newtonian fluids that exhibit both"viscous"and"elastic"characteristics in virtue of mechanisms to store energy and produce entropy. Usually the energy storage properties of such fluids are modelled using the same concepts as in the classical theory of nonlinear solids. Recently new models for elastic solids have been successfully developed by appealing to implicit constitutive relations, and these new models offer a different perspective to the old topic of elastic response of materials. In particular, a sub-class of implicit constitutive relations, namely relations wherein the left Cauchy-Green tensor is expressed as a function of stress is of interest. We show how to use this new perspective it the development of mathematical models for viscoelastic fluids, and we provide a discussion of the thermodynamic underpinnings of such models. We focus on the use of Gibbs free energy instead of the Helmholtz free energy, and using the standard Giesekus/Oldroyd-B models, we show how the alternative approach works in the case of well-known models. The proposed approach is straightforward to generalise to more complex setting wherein the classical approach might be impractical of even inapplicable.


Introduction
Most of the classical treatises on continuum mechanics when discussing constitutive theory start with the assumption that the stress is given as a function or a functional of the deformation, see for example Truesdell and Noll (2004), Müller (1985) and Coleman and Noll (1960). However, it has been argued that this approach is too restrictive, see especially Rajagopal (2003Rajagopal ( , 2006.
For example, it turns out that it is worthwhile to generalise the classical incompressible Navier-Stokes fluid model T = −pI + 2νD as (1.1a) where we use the standard notation T for the Cauchy stress tensor, p for the mechanical pressure, D for the symmetric part of the velocity gradient and f for a tensorial function. This approach to the generalisation of the classical incompressible Navier-Stokes fluid model is in contrast with the classical approach based on the formula T = −pI + g (D) , (1.2) where g is a tensorial function. (Here one can think of classical models for shear-thinning/thickening fluids.) Note that even if the relation D = f (T δ ) is invertible, and hence can be rewritten in the classical form T δ = g (D), it might be still convenient to use the representation (1.1b) because the representation (1.1b) might be much simpler than the classical one. (Meaning that the formula D = f (T δ ) might have a nice analytical form, while the inverse formula T δ = g (D) could be complicated or impossible to write down using elementary functions.) If (1.1b) is not invertible, then the two approaches are clearly strikingly different. Furthermore one can even go one step further and replace (1.1b) by a more general implicit relation k (T δ , D) = O, where k is a tensorial function and O denotes the zero tensor.
Similarly, in the context of mathematical models for elastic solids, one can generalise the standard isotropic Cauchy elastic materials wherein one assumes that T = h (B), where h is a tensorial function, and B denotes the left Cauchy-Green tensor. The generalisation of this formula is or even j (T, B) = O, where i and j are again tensorial functions. Clearly, everything that has been said in the context of models for fluids holds as well for the models for elastic solids.
Since viscoelastic fluids are non-Newtonian fluids that exhibit both "viscous" and "elastic" characteristics, it is desirable to see how to exploit the concept of implicit type constitutive relations in the context of mathematical models for these fluids. In particular, the phenomenological theory of viscoelastic rate-type fluids is based on the assumption that the elastic response of the fluid is modelled using the same concepts as in the classical theory of nonlinear solids. (Especially the assumed formulae for the Helmholtz free energy borrow heavily from the classical theory of nonlinear solids.) However, if one wants to model the elastic response in the spirit of recent advances in theory of nonlinear elastic solids that exploit constitutive relations of type (1.3), the phenomenological derivation of mathematical models viscoelastic rate-type fluids must be modified accordingly. Such modifications are discussed-with a special emphasis on the corresponding thermodynamic basis-in the current contribution.
The paper is organised as follows. First we briefly review the concept of constitutive relations for fluids and solids, that are developed using the approach (1.1) or (1.3) and its generalisations to the fully implicit setting k (T δ , D) = O or j (T, B) = O, see Section 2. Then we proceed with a brief review of the classical derivation of viscoelastic rate-type models that is based on the concept of natural configuration and the characterisation of the elastic response with the Helmholtz free energy, see Section 3.
In Section 4 we identify the Hencky strain tensor as the convenient tensorial quantity and the Gibbs free energy as the convenient thermodynamic potential that allows one to work with (elastic) constitutive response of type (1.3), and we outline a general procedure that allows one to develop novel viscoelastic rate-type models based on these concepts. In Section 5 we study the classical Giesekus/Oldroyd-B model both in the classical framework and the newly proposed framework. For this model it is possible to write down explicitly all formulae in both approaches, which allows us to clearly document the general theory. Finally, we also discuss the applicability of the maximisation of the entropy production hypothesis, see Rajagopal and Srinivasa (2004), in the theory of constitutive relations.

Materials specified via implicit constitutive relations
Before we proceed with viscoelastic fluids, we briefly summarise state-of-the-art regarding the novel approach to the mathematical modelling of viscous fluids and elastic solids.
2.1. Fluids. A simple generalisation of the Navier-Stokes fluid as indicated in (1.1) is the constitutive relation in the form D = α 1 + β T δ 2 n + δ T δ , (2.1) introduced in Le Roux and . (For a through discussion of a simpler model with δ = 0 see Málek et al. (2010).) With properly tuned parameter values this constitutive relation leads to the characteristic S-shaped curve in the strain rate-stress diagram, hence the model (2.1) can serve as simple phenomenological model for flows of substances that exhibit such behaviour, see for example Boltenhagen et al. (1997), Grob et al. (2014) and Mari et al. (2015) to name a few. (Further references and a thorough discussion can be found in Perlácová and Průša (2015) and Janečka and Průša (2015).) Note that if the constitutive relation leads to the S-shaped curve in the strain rate-stress diagram, then (2.1) is not invertible, and it can not be rewritten in the classical form T δ = g (D). Other models that are easy to describe using (1.1) are the models for incompressible fluids with pressure dependent viscosity, see Rajagopal (2006), or for that matter models wherein the viscosity depends on the stress such as the models used for the flow of ice, see Blatter (1995) and Pettit and Waddington (2003) to name a few and also  for further references regarding this class of models.
Other models that are easy to decscribe using (1.1) or the implicit relation k (T δ , D) = O are the models with activation criterion introduced by Bingham (1922) and Herschel and Bulkley (1926). These models are usually, see for example Duvaut and Lions (1976), written down in the form of a dichotomy relation where ν is a positive function and τ * is the yield stress.) If we use the alternative framework based on the constitutive relation k (T δ , D) = O, then these models can be rewritten as where x + = max{x, 0} and ν is a non-negative function. This observation can be exploited in the discussion of the physical and mathematical properties of the models, see for example Rajagopal and Srinivasa (2008) and Bulíček et al. (2012). For further discussion of generalised models of type k (T δ , D) = O we also refer the reader to Perlácová and Průša (2015). Note that in all cases discussed above the constitutive relation between the stress and the symmetric part of the velocity gradient has been related to entropy production mechanisms, see for example Janečka and Průša (2015) for a through discussion. Reader interested in analytical and semianalytical solutions to boundary value problems for fluids described by constitutive relations of the type T δ = g (D) or k (T δ , D) = O is referred to Málek et al. (2010), Le Roux and , Srinivasan and Karra (2015), Narayan and Rajagopal (2013), Fusi and Farina (2017), Fusi (2018), Housiadas et al. (2015), Gomez-Constante and  and Fetecau and Bridges (2020) to name a few. Numerical solution of the corresponding governing equations is investigated in Janečka et al. (2019), while rigorous numerical analysis for various models that fall into this class is discussed in Diening et al. (2013), Stebel (2016); Hirn et al. (2012), Tscherpel (2019), Farrell et al. (2020) and . Rigorous mathematical theory for some of the aforementioned models is developed in Bulíček et al. (2009Bulíček et al. ( , 2012 and Maringová andŽabenský (2018), see also Blechta et al. (2020) for a recent review.
Further generalisation of the algebraic relation k (T δ , D) = O is the implicit relation of the type that relates the histories of the stress tensor and the relative right Cauchy-Green tensor, for details see Průša and Rajagopal (2012). Special instances of (2.4) are rate-type and differential-type models for viscoelastic fluids and integral models for viscoelastic fluids, while in all these cases one can consider models with pressure/stress dependent material parameters, see for example Karra et al. (2011), Kannan and Rajagopal (2013), Housiadas (2015Housiadas ( , 2020 and Arcangioli et al. (2019) and references therein.

Solids.
Regarding the similar development in the case of constitutive relations for solids, wherein the classical approach is based on constitutive relations in the form T = h (B). The non-standard relations of the type B = i (T) or even j (T, B) = O have been considered without a proper thermodynamic basis or clearly articulated rationale for the choice of such constitutive relations in some early works such as Morgan (1966), Fitzgerald (1980) or Blume (1992) and Xiao and Chen (2003); Xiao et al. (2004). A systematic study of constitutive relations of the type B = i (T) with coherent explanation based on issues of causality for the choice of such constitutive relations was initiated by Rajagopal (2003), and the number of works focused on this concept has been growing ever since. Regarding elastic (non-dissipative) solids we refer the reader to a recent review by Bustamante and Rajagopal (2020) and references therein, while some aspects of the description of nonelastic response such as plasticity are discussed in Srinivasa (2016, 2015) and Cichra and Průša (2020). (A discussion and references on works dealing the invertibility of the constitutive relations can be found in Sfyris and Bustamante (2013), one of the earliest studies of the same is that by Truesdell and Moon (1975).) As an example of the benefits of the proposed change of perspective we can mention especially the study by Muliana et al. (2018), who address the classical topic of mathematical modelling of the mechanical response of rubber. From the current perspective, the important aspects of the novel approach to the constitutive theory of solids are the following. First, if one wants to work with constitutive relations of the type B = i (T), then it is convenient to work with the Hencky strain tensor H = def 1 2 ln B instead of the left Cauchy-Green tensor B. (This observation is made in various works in elasticity and inelasticity, see especially Xiao et al. (2004); Xiao and Chen (2003); Xiao et al. (2006). Note also that in the context of viscoelastic fluids the preferred quantity to work with is-from the numerical perspective-the log-conformation tensor, see Kupferman (2004, 2005); Hulsen et al. (2005), that is the logarithm of the conformation tensor. Apparently, this is not a coincidence.) Second, the thermodynamic potential of choice is the Gibbs free energy, see Rajagopal and Srinivasa (2011) and Srinivasa (2015). In particular, in the case of elastic solids the formula for the Hencky strain H = i (T) arises via the differentiation of a potential, and the potential can be identified as the Gibbs free energy, see especially Gokulnath et al. (2017) and Průša et al. (2020).

Viscoelastic rate-type fluids
Since viscoelastic fluids are fluids that exhibit both "viscous" and "elastic" characteristics, see for example Snoeijer et al. (2020) for a recent discussion concerning the same, one can ask whether the novel approach to the modelling of the response elastic solids can be also applied in this context. The point is the following. Traditionally the elastic properties of viscoelastic fluids are from the phenomenological perspective modelled using the same concepts as in the classical theory of nonlinear solids wherein the Cauchy stress tensor is assumed to be a function of the left Cauchy-Green tensor. The question is what happens if an opposite perspective to the classical one is taken-the left Cauchy-Green stress tensor is assumed to be a function the stress tensor, or one even assumes that these quantities are related by a implicit algebraic relation.
We focus on this question, and we show how to use the new perspective in the development of mathematical models for viscoelastic fluids, and discuss in detail the thermodynamic underpinnings of such models. In particular we focus on the use of Gibbs free energy instead of the Helmholtz free energy. Using the standard Giesekus/Oldroyd-B models, we show how the alternative approach works in the case of well-known models, and we argue that the proposed approach is straightforward to generalise to more complex settings.
Our basic framework will be the framework based on the concept of evolving natural configuration, see Eckart (1948) and Rajagopal and Srinivasa (2000), that has been successfully applied in various settings, see for example Málek et al. (2015b); Narayan et al. (2016); Málek et al. (2016), , Tůma et al. (2018), Málek et al. ( ),Řehoř et al. (2020 or Sumith et al. (2020) to name a few. (For other approaches to the mathematical modelling of viscoelastic rate-type fluids see Leonov (1976Leonov ( , 1987, Wapperom and Hulsen (1998) or the GENERIC framework, see Grmela andÖttinger (1997);Öttinger and  and Pavelka et al. (2018).) 4. Gibbs free energy based approach to viscoelastic rate-type fluids The concept of evolving natural configuration is based on the assumption that the overall response of a viscoelastic fluid is composed of a dissipative (viscous) response and a non-dissipative (elastic) response, while the latter can be described using the concepts known from the theory of nonlinear elasticity. However, the left Cauchy-Green tensor for the total response B must be replaced by the left Cauchy-Green tensor B κ p(t) associated to the elastic part of the fluid response. The subscript κ p(t) denotes the instantaneously relaxed or "natural" configuration, see Rajagopal and Srinivasa (2000) for a detailed rationale. (Another detailed explanation of the concept of natural configuration is given in a recent review article by . However, note that in the current contribution we do not explicitly use the decomposition of the total deformation gradient F = F κ p(t) G to the dissipative and non-dissipative part. In particular we do not explicitly work with the tensor G, which makes our approach different from the treatment in Rajagopal and Srinivasa (2000) and Málek et al. (2015a).) If we want to exploit the novel approach to the modelling of elastic response, we have to mimic (1.3), but the full Cauchy-Green tensor B must be replaced by the partial Cauchy-Green tensor B κ p(t) .
Furthermore, as we have already indicated, the partial left Cauchy-Green tensor B κ p(t) should be replaced by the corresponding Hencky strain tensor, and the thermodynamic potential of choice should be the Gibbs free energy. The details are worked out below. Note that unlike in other works focused on the use of Gibbs potential, see especially Narayan et al. (2015), we are not assuming that the processes of interest are isothermal. The theory outlined below is applicable in non-isothermal setting.
The main idea behind the presented approach is that the behaviour of the material in the processes of interest is determined by two factors, namely its ability to store energy and produce entropy. The energy storage mechanisms are in the present case specified by the choice of the Gibbs free energy g, while the entropy production mechanisms are specified by the choice of the formula for the entropy production ξ. Since the entropy η is obtained via differentiation of g, the evolution equation for η then necessarily links the assumed form of the Gibbs free energy g and the assumed entropy production ξ. This link can be exploited in the identification of the constitutive relations. Consequently, the first step of the proposed approach is to derive an evolution equation for the entropy. 4.1. Entropy evolution equation -general case. The fundamental equation in continuum mechanics is the evolution equation for the specific internal energy e, that is for the internal energy per unit mass, [e] = J kg reads see standard texts on continuum thermodynamics such as Truesdell and Noll (2004) or Müller (1985). (The notation in (4.1) is the standard one, the symbol T stands for the Cauchy stress tensor, D denotes the symmetric part of the velocity gradient, ρ denotes the density, and j q denotes the heat flux.) Evolution equation (4.1) is the starting point for the derivation of the sought entropy evolution equation.
In the present approach to the theory of viscoelastic rate-type fluids, we assume that the specific internal energy is a function of the specific entropy η, density ρ, and the Hencky strain tensor associated to the elastic response, Furthermore, it is convenient to split the Hencky strain tensor to the traceless (deviatoric) part H κ p(t) , δ and the spherical part h κ p(t) , where we introduce the notation A δ = def A − 1 3 (Tr A) I for the traceless (deviatoric) part of the corresponding tensor. We note that the deviatoric part of the Hencky strain tensor H κ p(t) , δ is equal to the rescaled Hencky strain tensor associated to the rescaled left Cauchy-Green tensor B κ p(t) that is defined as Průša et al. (2020), see also Xiao et al. (2004). The usage of the rescaled left Cauchy-Green tensor B κ p(t) is a common practice in the theory of slightly compressible solids, see for example Horgan and Saccomandi (2004). The rationale is that det B κ p(t) = 1, hence the deformation is conveniently split into a volume-preserving and volume-changing part.
We also note that if the elastic part of the response is volume-preserving, then det B κ p(t) = 1 and consequently Tr H κ p(t) = 0, which in this special case implies that Regarding the specific internal energy e we therefore assume that which allows us to seamlessly deal with incompressible materials as well. (If the material is homogeneous and incompressible, we can simply remove ρ from (4.6) since the density is in this case constant. Similarly, if the response of the elastic component is modelled as a response of an incompressible material, then h κ p(t) is a constant and we can simply remove h κ p(t) form (4.6). We however need to introduce the corresponding Lagrange multipliers that enforce the incompressibility constraint.) Introducing the thermodynamic temperature θ in the standard manner as we can rewrite (4.1) in terms of the Helmholtz free energy ψ, (In the next to the last term we can consider only the traceless part of ∂ψ ∂H κ p(t) , δ since it enters the dot product with the traceless tensor dH κ p(t) , δ dt .) Finally we can split the Cauchy stress tensor into the mean normal stress m = def 1 3 Tr T and the deviatoric part T δ , T = mI + T δ , (4.10) and we get We have also used the balance of mass dρ dt + ρ div v = 0.) We recall that the specific internal energy and the specific Helmholtz free energy are related by the Legendre transformation, and that the following formula holds (4.12) The manipulation outlined above therefore gives us the sought evolution equation (4.9) for the entropy η. This equation is the starting point for the specification of the constitutive relations. In the standard approach, one in principle wants to provide a constitutive relation for the Cauchy stress tensor T and formulae for the time derivatives of h κ p(t) and H κ p(t) , δ such that the right-hand side of (4.9) is nonnegative, which in turn guarantees that the second law of thermodynamics holds in the given material and during the given class of processes. (See for example Dostalík et al. (2019) for applications of this procedure to several well known viscoelastic rate-type models.) We however might not want to use (4.9) directly, we want to rewrite it in terms of Gibbs free energy g.
In order to use the Gibbs free energy we need to introduce the thermodynamic pressure p th , and the stress T κ p(t) , which we again split into the spherical part p κ p(t) and the traceless part T κ p(t) , δ , (4.13) Using the definitions we can further rewrite (4.9) in such a way that it is ready for the use of the specific Gibbs free energy g.
(The definition (4.14a) is the standard definition of the thermodynamic pressure as it is used in the classical thermodynamics of compressible fluids.) If we further introduce the notation the specific Gibbs free energy is then defined as We note that the transition to the Gibbs free energy is in fact done by another Legendre transformation, and that the use of Legendre transformation implies that (4.17d) If we use definitions (4.14), we see that (4.11) can be rewritten as which upon using formulae (4.17) yields This is a general evolution equation for the specific entropy η.

Entropy evolution equation -incompressible fluids.
Let us now manipulate (4.18) into a form convenient for the ongoing investigation of the constitutive relations. For the sake of simplicity we now restrict ourselves to incompressible fluids, which means that the density is a constant and it no longer enters the formula for the Gibbs free energy. Further we will not split the elastic response into the volume-preserving and volume-changing part, that is we use the full Hencky strain tensor H κ p(t) instead of the pair H κ p(t) , δ and h κ p(t) , and similarly for the reduced stress S κ p(t) . In this case the Gibbs free energy is therefore given as 20) and the evolution equation for the entropy (4.19) reduces to and relations (4.17) read We note that if we deal with an incompressible substance, then the Cauchy stress tensor is split to the spherical part −pI and the deviatoric part T δ , (4.23) and that the pressure p-which can be now understood as the Lagrange multiplier enforcing the incompressibility constraint-constitutes a new unknown field to be solved for. Namely, it is not given by an equation of state as a function of density, temperature and other variables. Furthermore the incompressibility constraint div v = 0 implies that D δ = D.
Finally, we note that if we deal with isotropic elastic response, which is the case in the rest of the paper, then (4.22b) in virtue of the standard representation theorem for isotropic functions, see for example Truesdell and Noll (2004), implies that the tensors H κ p(t) and S κ p(t) commute, that is 4.2.1. Manipulations with entropy production -direct use of Gibbs free energy. If we further assume that the Gibbs free energy has the property that the mixed derivative vanishes, that is if (4.25) 1 The notation might be slightly ambiguous here. Using the index notation we have where we have used the summation convention.
The reduced stress tensor S κ p(t) is a second order tensor in the current configuration that transforms normals to infinitesimal surfaces to traction vectors, hence we see that it is a tensor of type 2 0 . If we need to calculate the time derivative of such a tensor, then the natural concept is the upper convected derivative, (4.26) (See for example Stumpf and Hoppe (1997) for geometrical underpinnings of the concept of upper convected derivative and its link to the concept of Lie derivative.) Using the definition of the upper convected derivative, we see that (4.25) can be rewritten as 2 (4.27) (We are exploiting the fact that S κ p(t) commutes with ∂g ∂Sκ p(t) and the cyclic property of the trace. The fact follows again from the representation theorems for isotropic functions.) Finally, the standard manipulation of the heat flux term yields which is the formula used by Rajagopal and Srinivasa (2011).

4.2.2.
Manipulations with entropy production -indirect use of Gibbs free energy. Note however that (4.18) can be also manipulated differently, and in this case no further structural assumptions on the specific Gibbs free energy are necessary. In the incompressible setting (4.18) reads Unlike in the previous manipulation we keep the Hencky strain tensor H κ p(t) in the equation, but we replace its material time derivative by the upper convected derivative, (4.30) (The rationale for this manipulation is the same as in the case of the stress S κ p(t) .) Next, using the standard manipulation of the heat flux j q , we see that (4.18) can be finally rewritten as (4.31) 4.2.3. Constitutive relations. Either (4.31) or (4.28) can serve as a starting point for the specification of constitutive relations. Once we have postulated a formula for the Gibbs free energy, which is tantamount to the specification of the energy storage mechanism in the material of interest, we can proceed with the specification of the entropy production mechanisms. This boils down to the specification of the formula for the right-hand side of the generic entropy evolution equation that has the structure where j η denotes the entropy flux and ξ denotes the entropy production. Once the required entropy productionthe right-hand side or entropy evolution equation-is specified, we compare it with the entropy production implied by the choice of the Gibbs free energy, that is with the genuine right hand-side of (4.31) or (4.28), 2 For the sake of clarity we can resort to the index notation. We set where we use the summation convention. Since Sκ p(t) is a symmetric tensor that commutes with ∂g ∂Sκ , we see that Sκ p(t) is also a symmetric tensor.
which will allow us to fix the constitutive equations for T and T th . Examples of this procedure are given below in Section 5.2.
A more sophisticated version of this approach is based on the assumption of the maximisation of the entropy production, see Rajagopal and Srinivasa (2004). This approach is worked out in Section 5.3.

Temperature evolution equation.
Finally, let us remark that once the specific Gibbs free energy is given, one can then exploit the relation between the entropy and the specific Gibbs free energy, see (4.17a), and convert (4.28) or (4.31) respectively to an evolution equation for the temperature, indeed Note that if we use the assumption (4.24) that the mixed derivative of the specific Gibbs free energy vanishes, then the formula (4.33) simplifies to dη dt = − ∂ 2 g ∂θ 2 θ, S κ p(t) dθ dt . If it is not the case, we need to consider the term as well. The explicit formula for this term however depends on the rate-type constitutive relation for the stress tensor T κ p(t) .
5. Example -Giesekus/Oldroyd-B viscoelastic rate-type fluid in the approach based on the Gibbs free energy We now use the approach based on assuming the Gibbs free energy and appropriate entropy producing mechanisms to develop the popular incompressible rate-type viscoelastic fluid model due to Giesekus, see Giesekus (1982). (Note that this model reduces, for a special choice of parameters, to the Oldroyd-B model, see Oldroyd (1950).) In this case we have the luxury of writing down simple explicit formulas both for the Helmholtz free energy and the Gibbs free energy as well as for the entropy production. This makes the model ideal for clarifying the key ideas. 5.1. Giesekus/Oldroyd-B model via the Helmholtz free energy -classical approach. Let us first summarise basic facts regarding the classical incompressible viscoelastic rate-type model derived by Giesekus for heat conducting fluids with constant specific heat at constant volume c V,ref . (Note that the original derivation, see Giesekus (1982), is done in the purely mechanical context; no temperature evolution equation is given. The corresponding temperature evolution equation is derived for example in Wapperom and Hulsen (1998).) The governing equations for the mechanical quantities are in this case div v = 0, (5.1a) (5.1c) and the temperature evolution equation takes the form The Cauchy stress tensor T is given by the formulae where α ∈ [0, 1] is a model parameter, ν, ν 1 are positive material constants or functions that can eventually depend on temperature, and µ is a material constant or a function proportional to the thermodynamic temperature; µ plays the role of the shear modulus for the elastic response, while the combination ν1 µ(θ) plays the role of relaxation time, see (5.1c). (More complex temperature dependence of µ is not allowed. If we were dealing a more complex function than the constant/linear function of temperature, then we would get a model with specific heat at constant volume that depends on B κ p(t) , see for example comments in Hron et al. (2017). While such models might be useful, we do not consider them in the present contribution.) As we shall see later the material constant/function µ enters, unlike ν and ν 1 , the formula for the specific Helmholtz free energy ψ, whose derivatives with respect to θ and B κ p(t) determine the properties of the material. Therefore we shall explicitly write µ(θ) in order to emphasise the temperature dependence of µ. On the other hand, for ν and ν 1 the temperature dependence is not-from the theoretical point of view-as important as in the case of µ, hence we shall simply write ν and ν 1 instead of ν(θ) and ν 1 (θ). Finally, the symbol b denotes the body force.
A few remarks concerning the system (5.1) are in order. First, we note that the right-hand side of (5.1c) can be rewritten as − µ(θ) αB 2 In this form it is straightforward to see that if B κ p(t) = I, then the right-hand side of (5.1c) vanishes. (This is a useful information if one is interested in steady states predicted by (5.1) in the case of zero external forcing.) Second, the last term in the temperature evolution equation (5.1d) can be rewritten as This manipulation shows that the corresponding term is for α ∈ [0, 1] nonnegative, and that it vanishes for B κ p(t) = I, which a prospective steady state in the system without external forcing. Third, we recall that if we set α = 0, then we obtain the standard Oldroyd-B model, see Oldroyd (1950). The specific Helmholtz free energy ψ is in this case of model (5.1) known to be given by the formula and the entropy production reads Note that (5.3) shows that the entropy production is nonnegative. We recall that the Hencky strain tensor associated to the elastic response H κ p(t) is given by (4.2), hence we have B κ p(t) = e 2Hκ p(t) , (5.6) and the formula for the Helmholtz free energy ψ in terms of H κ p(t) reads ψ(θ, H κ p(t) ) = def ψ thermal (θ) + µ(θ) 2ρ Tr e 2Hκ p(t) − 3 − ln det e 2Hκ p(t) . (5.7) Using the fact that the Hencky strain tensor is represented by a symmetric matrix and the standard identity det e A = e Tr A , we see that (5.7) can be rewritten as Tr e 2Hκ p(t) − 3 − 2 Tr H κ p(t) . (5.8) The stress T κ p(t) is then given via the definition , see (4.14), which in our case yields where we have used Lemma 1. Equation (5.9) can be solved explicitly for the Hencky strain H κ p(t) , Note that if we rewrite (5.9) in terms of B κ p(t) , then we get Having obtained formula (5.10) we can use the definition of the specific Gibbs free energy g = ψ − and we can write down an explicit formula for the Gibbs free energy corresponding to the Helmholtz free energy introduced in (5.4), which is easy to convert to the final form Since the Legendre transformation from the Helmholtz free energy to the Gibbs free energy is done with respect to the deformation and the temperature variable is kept intact, we see that the function g thermal (θ) is identical to the function ψ thermal (θ).
We note that if we use spectral representation of S κ p(t) , then (5.12) can be rewritten as . The function f (x) = x − (x + 1) ln(x + 1) that appears on the right-hand side of (5.14) is a nonpositive concave function that vanishes if and only if x = 0. Furthermore, by appealing to (5.13) we see that the Gibbs free energy becomes infinite once the quantity Tκ p(t) µ(θ) + I vanishes. This means that the quantity Tκ p(t) µ(θ) + I remains a positive definite matrix. (This is equivalent to the positive definiteness of B κ p(t) in the classical approach.) Finally, we observe that the Gibbs free energy vanishes for T κ p(t) = O, which corresponds to the fact that the Helholtz free energy vanishes for B κ p(t) = I.

Approach based on Gibbs free energy.
Let us now assume that the Gibbs free energy g is given by (5.13), and let us rederive the governing equations (5.1) within the approach based on the Gibbs free energy. (Instead of the Helmholtz free energy we start with the Gibbs free energy, and we observe what steps need to be taken, if we want to rederive the Giesekus model.) Using the standard matrix calculus we can observe that which is the expected result, since the use of Legendre transformation implies that H κ p(t) = − ∂g ∂Sκ p(t) (θ, S κ p(t) ).
Compare (5.15) and (5.10). Now we need to exploit the evolution equation for the entropy. In particular, we need to propose an evolution equation for T κ p(t) and a constitutive relation for T in such a way that the entropy production is nonnegative. (The thermal part-the heat flux-can be manipulated in the standard manner.) can be expressed as while the motivation for this manipulation is to get the term T κ p(t) • D, which can be conveniently manipulated in (5.16). In (5.17) we are using the cyclic property of the trace and the fact that T κ p(t) and H κ p(t) commute, which implies that T κ p(t) and e −2Hκ p(t) commute as well. If we substitute (5.17) into (5.16), we get The primitive quantity in the approach based on the Gibbs free energy is the reduced stress tensor S κ p(t) = Tκ p(t) ρ , and we need to find an evolution equation for this quantity. (Or for that matters an evolution equation for T κ p(t) since these quantities differ only by the multiplication by the density ρ, which is in our case a constant.) We recall that the evolution equation for T κ p(t) and the constitutive relation for T must be chosen in such a way that the corresponding terms on the right-hand side of (5.16) are nonnegative. 5.2.1. Oldroyd-B model. Before we discuss the Giesekus model, we focus on the case α = 0 which leads to the Oldroyd-B model. Regarding the entropy production terms in (5.18) the following simple choice (5.19b) is of particular interest. If we use (5.19) in (5.18), then we get (5.20) hence the entropy production in mechanical processes in nonnegative as required. Indeed, the term in the entropy production can be rewritten as , where we have used the fact that H κ p(t) = − ∂g ∂Sκ p(t) (θ, S κ p(t) ) and the particular formula for the Gibbs free energy, see (5.15). Note that the last observation implies that if we use constitutive assumptions (5.19) and (5.15), then the entropy evolution equation (5.20) rewritten in terms of T κ p(t) reads (5.23) In this case the nonnegativity of the entropy production in mechanical processes still holds-equation (5.23) is the same as equation (5.20) but rewritten in terms of different variables. However, the nonnegativity of the entropy production is less straightforward to recognise than in the formulation based on H κ p(t) .
If we further assume that the heat flux is given by the standard Fourier law, where κ is a nonnegative constant, then the entropy production both in thermal and mechanical processes is nonnegative and the second law of thermodynamics holds. Once we have identified the energy storage mechanisms-specific Gibbs free energy (5.13)-and the entropy production mechanisms-terms on the right-hand side of the entropy evolution equation (5.20)-the complete system of governing equations is obtained as follows. First, we use (5.19a) and (5.9), and we can conclude that the evolution equation for H κ p(t) reads (5.25) which can be as well rewritten in terms of T κ p(t) as (5.26) (In the latter case it suffices to take the upper convected derivative of (5.9) and use (5.19a). We recall that the definition of the upper convected derivative implies that ▽ I = −2D.) Finally, if we wanted to, we could go back to the left Cauchy-Green tensor B κ p(t) , see (5.6), which allows us to rewrite (5.25) as (5.27) This is the evolution equation (5.1c) for α = 0. If we use (5.26), then the system of governing equations for the mechanical quantities reads div v = 0, (5.28a) where the Cauchy stress tensor T is given by the formula Since (5.28c) can be rewritten as (5.27), and since (5.9) holds, that is T κ p(t) = µ(θ) e 2Hκ p(t) − I holds, we see that (5.28) is equivalent to the mechanical part of (5.1). (A minor redefinition of the mean normal stress m is necessary. This is however an inconsequential modification of the governing equations.) Regarding the temperature evolution equation (5.1d) we have to start with the entropy evolution equation either in the form (5.20) or in the form (5.23). In particular, we need to show that the entropy production identified in (5.20) or equivalently in (5.23) is tantamount to the entropy production (5.5) that leads to the Oldroyd-B model. (We recall that Oldroyd-B model corresponds to α = 0.) The critical term can be manipulated as follows Making use of identity (5.3) we see that the entropy production mechanisms postulated in (5.20) or (5.23) are identical to those postulated in the standard derivation of the Oldroyd-B model (5.5). Note that (5.29) explicitly shows that the entropy production for the Oldroyd-B fluid can be equivalently written down either in terms of a algebraic formula for B κ p(t) or in terms of algebraic formula to the rates ▽ B κ p(t) . Since the entropy production mechanisms ξ are now confirmed to be the same in both approaches, we see that the right-hand side of the entropy evolution equation is the same both for the model derived via the Gibbs free energy and the Helmholtz free energy. Finally, it suffices to recall that the the entropy η is given in terms of derivatives of the potentials as η(θ, S κ p(t) ) = − ∂g ∂θ (θ, S κ p(t) ) and η(θ, B κ p(t) ) = − ∂ψ ∂θ (θ, B κ p(t) ), see (4.22a) and (4.12). Consequently, if we substitute for η into (5.30), we get for both approaches the same evolution equation for the temperature. (In the case of the Gibbs free energy the primitive variable will be the reduced stress tensor S κ p(t) , while in the case of Helmholtz free energy the primitive variable will be the left Cauchy-Green tensor associated to the elastic response B κ p(t) .) For readers convenience both approaches are compared in Summary 1 and Summary 2. 5.2.2. Giesekus model. Regarding the Giesekus model, the derivation in the approach based on the Gibbs free energy is a slight modification of the derivation of the Oldroyd-B model presented above. The Gibbs free energy remains the same, see (5.13), and the entropy evolution equation remains the same (5.18) as well, that is The entropy production mechanisms are however modified to Reduced stress (notation): ρ Entropy production: Material parameters: c V,ref specific heat at constant volume -positive constant; µ(θ) shear modulusnonnegative function, typically proportional to θ; ν, ν 1 viscosity -nonnegative functions of the primitive variables, typically constants; κ thermal conductivity -nonnegative function of the primitive variables, typically constant Evolution equations (mechanical variables m, v, T κ p(t) and thermal variable θ): Cauchy stress tensor: T = mI + 2νD δ + T κ p(t) , δ Thermodynamical relations: Left Cauchy-Green tensor B κ p(t) and Hencky strain tensor H κ p(t) associated to the elastic response: which implies constitutive relations

Summary 2: Incompressible Oldroyd-B model via Helmholtz free energy
Specific Helmholtz free energy: Entropy production: Material parameters: c V,ref specific heat at constant volume -positive constant; µ(θ) shear modulusnonnegative function, typically proportional to θ; ν, ν 1 viscosity -nonnegative functions of the primitive variables, typically constants; κ thermal conductivity -nonnegative function of the primitive variables, typically constant Evolution equations (mechanical variables m, v, B κ p(t) and thermal variable θ): Cauchy stress tensor: Left Cauchy-Green tensor B κ p(t) and Hencky strain tensor H κ p(t) associated to the elastic response: It is straightforward to check that the entropy production specified in (5.32) is nonnegative. Indeed, the critical term can be manipulated as follows where we have used the relation (5.10) between the stress tensor T κ p(t) and the Hencky strain tensor H κ p(t) . The first term on the right-hand side of (5.34) is up to a positive factor (1 − α) the same as the term in the entropy production for the see (5.23), while the nonnegativity of the term αT κ p(t) • T κ p(t) is evident. (Note that the entropy production is particularly simple if we set α = 1.) Using the same steps as in the previous section it is then straightforward to check that the constitutive relations (5.33) lead, once rewritten in terms of B κ p(t) , to the governing equations for the standard Giesekus viscoelastic rate-type fluid.
To conclude, we see that the Giesekus viscoelastic rate-type fluid is specified by the Gibbs free energy in the form ρ , while the entropy production takes the form (5.37) 5.3. Approach based on Gibbs free energy and entropy production maximisation. In a nutshell the hypothesis of entropy production maximisation, see Rajagopal and Srinivasa (2004), states that once the entropy production ability of the material is specified in terms of entropy production ξ given as a function of thermodynamics affinities/fluxes, and once all other structural constraints are taken into account, then the constitutive relations for fluxes/affinities as functions of affinites/fluxes must be chosen in such a way that the entropy production is by this choice of fluxes/affinities maximised. (See also  for an explanatory presentation and worked out examples.) In the remainder of this section we show how the procedure works for the Giesekus/Oldroyd-B models in the approach based on the Gibbs free energy. Note however that the terminology affinity/flux, albeit the standard one, might be seen as inappropriate for several reasons, see Rajagopal and Srinivasa (2019). In the current contribution we therefore take the liberty to use the term resultant instead of affinity and the term determinant instead of flux, see Rajagopal and Srinivasa (2019) for the rationale of this terminology. 5.3.1. Oldroyd-B model. The Gibbs free energy remains the same as in the preceding sections, see (5.13), and the entropy evolution equation remains the same (5.18) as well, that is In virtue of thermodynamic relation H κ p(t) = − ∂g ∂Sκ p(t) (θ, S κ p(t) ) we for the given Gibbs free energy (5.13) get T κ p(t) = µ(θ) e 2Hκ p(t) − I , (5.39) see Section 5.2 for details. This relation allows us to rewrite (5.38) as (5.40) This is the entropy production as it is implied by the choice of the Gibbs free energy (5.13). (Note that (5.40) is unlike (5.38) written down in terms of the stress tensor T κ p(t) .) The rate-type quantities and can be identified as resultants, while the corresponding quantities in the dot products in (5.40) can be identified as determinants, (5.42b) (We also have the classical determinant-resultant pair j q and ∇θ.) On the other hand the assumption regarding the specification of the entropy production is for the Oldroyd-B model the following (5.43) see (5.23) and the subsequent discussion. (We need to rewrite (5.37) in such a way that it contains the resultants identified in (5.41).) The constitutive relations we need to identify are the constitutive relations for the Cauchy stress tensor T and the stress tensor T κ p(t) and the heat flux j q . The entropy production maximisation procedure requires us to maximise the entropy production ξ given by (5.43) with respect to resultants (5.41) and ∇θ, and subject to the constraint (5.44) that is implied by the right-hand side of (5.40). The constrained maximisation problem therefore reads where we have used the determinant/resultant notation introduced in (5.41) and (5.42). The auxiliary functional for the constrained maximisation problem reads Let us first solve (5.47) for the Lagrange multiplier λ. Once we identify the value of the Lagrange multiplier λ, it will be straightforward to find the determinants j q , J 1 and J 2 for which the entropy production attains its maximum. Taking the dot product of individual equations in (5.47) with the corresponding resultant, and taking the sum of these equations yields Making use of constraint (5.44), we see that the right-hand side of (5.48) is equal to ξ. Furthermore, the direct differentiation of ξ given by the formula (5.43) yields (This is not surprising since the entropy production is "quadratic" in resultants.) Making use of (5.50) in (5.48) allows us to easily identify the Lagrange multiplier λ. Indeed (5.48) implies that 2 1+λ λ ξ = ξ, hence 1 + λ λ = 1 2 . (5.51) Having obtained the formula for the Lagrange multiplier, we can go back to (5.47) and find formulas for the determinants. We see that equations (5.47) immediately imply that (5.52c) If (5.52c) holds, then the matrices and ▽ Tκ p(t) µ(θ) + I commute, which means that (5.52c) in fact simplifies to This is the evolution equation for T κ p(t) in the case of Oldroyd-B model, see (5.26) and the discussion following this equation.
In order to show that the matrices µ(θ) + I associated to the eigenvalueλ. (Since the matrix is symmetric positive definite we know that the normalised eigenvectors form an orthonormal basis, and that the eigenvalues are positive. This also guarantees that the matrix is invertible.) Such an eigenvector is also an eigenvector of matrix T κ p(t) with eigenvalue µ(θ) λ − 1 . Making use of (5.52c) we see that The only difference is in the postulated entropy production mechanisms. Instead of entropy production ξ given by (5.43) we now set (5.56) Few observations regarding (5.56) are at hand. First, if we set α = 0, then we recover the postulated entropy production ξ for the see (5.43). Second, since the matrix Tκ p(t) µ(θ) + I is positive definite, and since α ∈ [0, 1], we see that the postulated entropy production is nonnegative. In fact, the critical term in the entropy production has the structure of a weighted dot product of matrices ▽ Tκ p(t) µ(θ) + I that represent the resultant A 2 , see (5.41b).
The constrained maximisation problem (5.45) is now solved by the same procedure as in the case of Oldroyd-B model. (The identification of determinants and resultants remains the same.) In particular, the "quadratic" structure of the entropy production allows one to easily identify the Lagrange multiplier λ as a solution to the equation 1+λ λ = 1 2 . The counterpart of (5.52c) is then (5.57) and regarding the matrices on the left-hand side, it can be shown again that they commute. Equation (5.57) therefore collapses to which can be rewritten as (5.60) Simple substitution of (5.11) then reveals that (5.60) is tantamount to the standard evolution equation for B κ p(t) used in the Giesekus model, see (5.1c) and (5.2).

Conclusion
We have shown how a novel approach to the constitutive relations for elastic solids can be embedded in the development of mathematical models for viscoelastic fluids, and we have discussed in detail thermodynamic underpinning of such models. In particular, we have focused on the use of Gibbs free energy instead of the Helmholtz free energy. Using the standard Giesekus/Oldroyd-B models as examples, we have show how the alternative approach works in the case of well-known models and in the fully non-isothermal setting. The proposed approach is straightforward to generalise to more complex setting wherein the classical approach based on the Gibbs potential might be impractical of even inapplicable.

Appendix A. Daleckii-Krein formula and its consequences
Since we frequently need to differentiate matrix valued functions, it is worth recalling that the derivatives of matrix valued functions are easy to obtain using the Daleckii-Krein formula, see Daletskii and Krein (1965). (Short proofs of the same are given in Bhatia (2015, Theorem 5.3.1) and Bhatia (2013, Theorem V.3.3).) Regarding the matrix valued functions and especially tools for computations thereof we refer the reader to Higham (2008).
Theorem 1 (Daleckii-Krein formula). Let A is a real symmetric matrix in R k×k with spectral decomposition are the eigenvalues of A, and {P i } k i=1 denotes the projection operator to the i-th (normalised) eigenvector v i , that is P i = def v i ⊗ v i . Let f is a continuously differentiable real function defined on the open set containing the spectrum of A. Then the corresponding matrix valued function f, is differentiable, and the Gateaux derivative of f at point A in the direction X is given by the formula Furthermore, let [A ○ B] ij = def A ij B ij denotes the elementwise Schur/Hadamard product of matrices A and B.
(No summation convention.) Using the Schur product we can rewrite (A.2) as provided that the Schur product is taken in a basis in which A is diagonal, and for i = j the quotient in (A.3) is interpreted as f ′ (λ i ).
Note that in the continuum mechanics literature the same formulae have been in special cases rediscovered independently, see form example Jog (2008) for the case of matrix logarithm. Using the Daleckii-Krein formula it is straightforward to show the following two lemmas we have been frequently using in our analysis. where f ′ is the matrix valued function associated to f ′ , that is where A = ∑ k i=1 λ i P i denotes the spectral decomposition of A.
Proof. Since the matrices A and B are symmetric matrices and they commute, they share the same set of eigenvectors {v i } Lemma 2. Suppose that A and B are real symmetric matrices in R k×k that are differentiable functions with respect to t ∈ R. Furthermore, suppose that these matrices commute for all t ∈ R, that is AB = BA, and let f be a continuously differentiable real function defined on the open set containing the spectrum of A at every t ∈ R. Then where f ′ is the matrix valued function associated to f ′ , see (A.5).
Proof. We use the chain rule, Daleckii-Krein formula and spectral decompostion A = ∑ k i=1 λ i P i and B = ∑ k i=1 µ i P j , where we have used the fact that P i P j = δ ij P j and the properties of trace. (No summation.) Note that the identity does not in general hold without the trace. For example it is not in general true that de A dt = e A dA dt .