A Simple Construction of a Thermodynamically Consistent Mathematical Model for Non-Isothermal Flows of Dilute Compressible Polymeric Fluids

: We revisit some classical models for dilute polymeric ﬂuids, and we show that thermodynamically consistent models for non-isothermal ﬂows of these ﬂuids can be derived in a very elementary manner. Our approach is based on the identiﬁcation of energy storage mechanisms and entropy production mechanisms in the ﬂuid of interest, which, in turn, leads to explicit formulae for the Cauchy stress tensor and for all of the ﬂuxes involved. Having identiﬁed these mechanisms and derived the governing equations, we document the potential use of the thermodynamic basis of the model in a rudimentary stability analysis. In particular, we focus on ﬁnite amplitude (nonlinear) stability of a stationary spatially homogeneous state in a thermodynamically isolated system.


Introduction
Starting from the seminal work by Kramers [1] kinetic-type models have been widely used in the mathematical modelling of polymeric fluids, see the monographs by Bird et al. [2], Beris and Edwards [3], Öttinger [4], Huilgol and Phan-Thien [5], Dressler et al. [6], Öttinger [7] Kröger [8], and the review paper by Lozinski et al. [9] to name a few (a mathematically inclined reader-a specialist in theory of partial differential equations-is also referred to the concise presentation in Le Bris and Lelièvre [10]). The overwhelming majority of the works, especially in the field of numerical simulations, see, for example Lozinski and Chauvière [11], Knezevic and Süli [12] or Mizerová and She [13], are however restricted to isothermal flows. In particular, the temperature field is often tacitly assumed to be homogeneous in space, and an evolution equation describing the temporal and spatial variations of the temperature field is rarely formulated, albeit some approaches, such as the GENERIC formalism, see Öttinger [7], or Pavelka et al. [14], allow one to do so. (regarding the classical macroscopic models, temperature evolution equations have been formulated for example in Peters and Baaijens [15], Wapperom and Hulsen [16] and Dressler et al. [6], see also Hron et al. [17] for further discussion). In the present contribution, we provide a straightforward self-contained derivation of a simple kinetic-type model for non-isothermal flows of compressible dilute polymeric fluids. The model is essentially the same as the model discussed in Öttinger and Grmela [18], but the approach to the derivation of a thermodynamically consistent model is substantially different.
We only use elementary arguments easily accessible to everyone familiar with basic principles of continuum mechanics such as the general form of balance equations. It turns out that, once one correctly specifies the energy storage mechanisms and the entropy production mechanisms in the fluid of interest, then some formulae, in particular the Kramers formula for the stress tensor, are natural consequences of the basic principles (there is no need to assume these formulae a priori). The same holds true in the derivation proposed by Öttinger and Grmela [18], where the authors claim that this fact is a "quite remarkable" consequence of the GENERIC structure. While the GENERIC framework has definitely its merits, we show that much more elementary arguments lead-regarding this issue-to the same conclusion. Furthermore, having obtained a thermodynamically consistent model, we exploit thermodynamic considerations in a rudimentary stability analysis in the spirit of Coleman [19], see also Bulíček et al. [20] for a discussion.
The paper is organised, as follows. First we recall basic principles of continuum mechanics, see Section 2, with particular attention placed on the basic tenets of the classical kinetic-type theory for dilute polymeric fluids. The key assumption regarding dilute polymeric fluids is that the polymeric chains do not effectively contribute to the mass of the solvent/polymer mixture, which means that the density of the mixture coincides with the density of the solvent. Furthermore, the kinetic energy of the polymeric chains is also neglected, the only mechanical energy contribution of the polymeric chains being via the energy stored in the stretched polymeric chains.
In Section 3, we focus on the dynamics of polymeric chains. Regarding the chains, we do not model them individually, but we again follow a kinetic-type theory, and we instead formulate a general form of the Fokker-Planck equation for the configurational distribution function of the polymeric chains. In our setting, the Fokker-Planck equation contains two unknown flux terms, namely the flux in configurational space and the flux in physical space. The latter is essential if we want to model the so-called centre-of-mass diffusion. Furthermore, we also comment on possible stationary solutions of the Fokker-Planck equation, and we carefully discuss the interplay between boundary conditions in the configurational space and asymptotic behaviour of the spring potential.
In Section 4, we give a formula for the specific Helmholtz free energy for the given fluid, which is we specify possible energy storage mechanisms. Subsequently, we proceed in Section 5 with the derivation of the constitutive relations, that is with the derivation of the equations relating the Cauchy stress tensor, the energy flux, the fluxes in the Fokker-Planck equation, and the kinematical quantities. The basic idea is that we look for closure relations that guarantee the nonnegativity of the entropy production. In this sense, the identification of the entropy production mechanisms and energy storage mechanisms provides a complete characterisation of the fluid of interest. A summary of the resulting model is given in Section 5.4.
Once we obtain the constitutive relations we discuss, see Section 6, the implications of the thermodynamic basis regarding the investigation of the stability of container flows (from the thermodynamic point of view, the flow in a thermally and mechanically isolated container constitutes an example of a thermodynamically isolated system). We show that the system of governing equations has a "trivial" stationary solution, and using thermodynamic arguments, we explicitly construct a nonnegative functional that decays in time and vanishes if and only if the system reaches the spatially homogeneous stationary state. The construction of such a functional is clearly a precursor for a rigorous stability analysis, which is, however, beyond the scope of the present contribution.

Preliminaries
The classical balance equations for the mass and momentum and the evolution equation for the internal energy of a continuous medium are The symbol ρ denotes the density of the medium, v denotes the Eulerian velocity field, ⊺ denotes the symmetric velocity gradient, T denotes the Cauchy stress tensor, e denotes the specific internal energy, j e denotes the energy flux, and b denotes the external body force.
denotes the material derivative, and the symbol A ∶ B = def Tr (AB ⊺ ) denotes the Frobenius norm.
The subscript x reminds us that the differential operator concerned is applied with respect to the spatial variable x. Regarding the derivation of the balance equations, see, for example, Málek and Průša [21] or any standard textbook on continuum mechanics, such as Müller [22] or Gurtin et al. [23]. The basic idea in a kinetic-type theory of dilute polymeric fluids is that the mass contribution of the polymer chains is assumed to be negligible, hence we can replace the total density ρ in (1) by the solvent density ρ s . Moreover, the barycentric velocity of the solvent/polymer mixture coincides with the velocity of the solvent, hence v in (1) can be interpreted as the velocity of the solvent only. The only place where the polymeric chains enter the balance equations is the formula for the internal energy e, the Cauchy stress tensor T and the energy flux j e . All of these assumptions/simplifications are common in kinetic-type theories of dilute polymeric fluids, and we also adopt them in the current contribution.
For further reference, we can therefore write the balance equations as and we recall that the unknown fields ρ s and v are functions of the spatial position x and time t. We note that, if one wanted to develop a model wherein the mass of the polymeric chains is not negligible, then one would need to use a variant of mixture theory, see, for example, Rajagopal and Tao [24] or Hutter and Jöhnk [25], and also Souček et al. [26] for a careful discussion.

Fokker-Planck Equation
We are in a position to describe the dynamics of the polymeric chains dispersed in the solvent. In what follows, we deal with a very simple setting, and we assume that the whole polymeric chain can be modelled as a single dumbbell (two beads connected with a spring), while more sophisticated models are typically necessary in order to obtain quantitative agreement between the model predictions and the experimental data (see the aforementioned monographs for a list of such models. Additional references regarding viscoelastic models can be also found in Vinogradov and Malkin [27], Larson [28] or Leonov and Prokunin [29]. A historical perspective is given in Tanner and Walters [30]). Note however that a generalisation to more complex polymeric chains/more sophisticated models is straightforward.
The dynamics of the individual polymeric chain are governed by a stochastic differential equation, which has the associated Fokker-Planck equation, see, for example, Le Bris and Lelièvre [10] or Barrett and Süli [31] for details. If we consider dumbbells, their configuration is fully determined by the end-to-end vector q, and the Fokker-Planck equation governs the evolution of the configurational distribution function ϕ(t, x, q) in the configurational space D of all admissible end-to-end vectors. The physical meaning of the configurational distribution function is clear, the integral gives the number of polymer chains at time t and spatial location x whose configuration vector lies in the domain D ′ . The integral over the whole D gives the particle number density n p (t, x), and the integral over a spatial domain gives one the number of polymeric chains in the spatial domain Ω ′ .
In the case of Hookean dumbbells, the configurational space D is the whole R 3 , since the stretch of the spring is not limited. In the case of the FENE dumbbell model (finitely extensible nonlinear elastic spring) the configurational space D is a ball in R 3 centred at the origin with a given radius, and the radius determines the maximum stretch of the spring.

Fokker-Planck Equation in the Case of Velocity Field with Nonzero Divergence
The Fokker-Planck equation has, in our case, the form where j ϕ,x and j ϕ,q are unknown flux terms. Our task will be to identify these flux terms using the known energy storage mechanisms and entropy production mechanisms. We note that the spatially dependent term reads div x (vϕ) and not vdiv x ϕ (this subtle difference does not matter in the case of incompressible fluids, but it is critical provided that we are working with a compressible fluid. See. for example, Degond and Liu [32] for a careful discussion of the Fokker-Planck equation). For further reference we also record that (7) can be rewritten as

Boundary Condition in the Configurational Space
The boundary condition on ∂D is the no-flux boundary condition where n q denotes the unit outward normal to the set of admissible end-to-end vectors D ⊂ R 3 (in the case of Hookean dumbbells where D = R 3 the boundary condition is interpreted as a decay of the corresponding function at infinity). This choice of boundary condition has several implications regarding the evolution equation for polymer number density, see below, and regarding the behaviour of the configurational distribution function on the boundary of D, see Section 3.4.

Evolution Equation for Polymer Number Density
The evolution of the polymer number density n p introduced in (5) is governed by the partial differential equation This equation is straightforward to obtain via integration of the Fokker-Plank equation with respect to the configurational variable q. The boundary term in the configurational space vanishes by virtue of (9). The last equation can be rewritten as where d dt denotes the standard material time derivative. Equation (11) shows that if j ϕ,x = 0 and if the velocity field does not identically satisfy div x v = 0, then, in general, we cannot expect that the polymer number density n p at the given material point is preserved. This makes the study of the dynamics of a compressible fluid with polymeric centre-of-mass diffusion substantially different from the simpler case of an incompressible fluid without polymeric centre-of-mass diffusion, see, for example, Barrett and Süli [33] and Feireisl et al. [34].

Force Potential
The force F in the spring connecting the beads is determined by a spherically symmetric potential where q ref is a characteristic length scale for the spring. (Note the sign convention in (12) and the sign convention used later in (15).) For example, a popular choice for the force potential is the potential that is introduced by Warner [35] for the FENE dumbbell model. The potential is given by the formula We note that the choice (13) predicts infinite force as s → b 2 −. This means that an infinite force is necessary to maximally stretch the spring, which is the expected behaviour in view of the fact that we are dealing with a finitely extensible spring. The spring can not be extended beyond a certain length no matter what force is applied.
In what follows, we will not work with a specific force potential, we will only assume that the force potential U is composed of two parts U = U e + U η,θ , where and where θ ref is a reference temperature, and q ref is a characteristic length scale. The potential U η,θ is the "entropic" part of the spring potential, and it is proportional to the ambient temperature θ s , while the potential U e is independent of the temperature (see, for example, Ericksen [36] for the rationale of the nomenclature, and the fact that it makes sense to consider the potential proportional to the temperature). If we are dealing with a finitely extensible spring, then the force tends to infinity as the spring reaches its maximum length (see for example the FENE potential introduced in (13)).
This suggests that the value of the configurational distribution function ϕ on the boundary of D should be zero. We informally show that this property is in fact encoded in the no-flux boundary condition (9). If we for a moment suppose that the Fokker-Planck Equation (7) reads that is if we consider specific formulae for the fluxes j ϕ,q and j ϕ,x , then (9) reduces to Here, θ s denotes the solvent temperature, ζ denotes the hydrodynamic drag coefficient and k B denotes the Boltzmann constant. If D is a ball in R 3 , then q ∂D = q max n q , and (16) reduces to Furthermore, in the case of spherically symmetric potential, the force F is parallel to q, and if we are dealing with a finitely extensible spring, then the potential is chosen in such a way that the magnitude of F approaches infinity as the vector q approaches the boundary of D. Having this piece of information in hand we inspect (17), and we see that if (17) holds, then the middle term 2F ζ ϕ • n q constrains the configurational distribution function ϕ to decay to zero faster than the growth of the force F as one approaches the boundary of D. Consequently, we can claim that (17) implies ϕ ∂D = 0.

Stationary Solution of the Fokker-Planck Equation in a Spatially Homogeneous State at a Given Temperature
If we assume for a moment that the fluxes in the general Fokker-Planck equation have been successfully identified, and have the form discussed in (15), then we can explicitly identify the stationary solution to the Fokker-Planck equation in the case of a spatially homogeneous temperature field and configurational distribution function field (partial derivatives with respect to time and the spatial variable x vanish). Under these assumptions, we see that the Fokker-Planck Equation (15) reduces to div q 2F If we assume that the force in the spring is given in terms of a potential, see (12), and if we substitute (12) into (18), then we get Using the formulae for the differential operators in a spherical coordinate system, we see that (19) reduces to a single ordinary differential equation where we have used the notation q = def q . We immediately see that the solution to (20) reads where C is a constant. The constant is fixed by the condition n p (t, x) = ∫ D ϕ(t, x, q) dq, with ϕ now independent of t and x here, and hence the configurational distribution function that solves the stationary Fokker-Planck in the case of spatially homogeneous fields is given by the formula ϕ = M n p ,θ s , where and n p is a positive constant. This configurational distribution function is the expected "equilibrium" configurational distribution function in the case when the dilute polymeric fluid occupies a thermodynamically isolated vessel. We return to this issue in Section 6.

Helmholtz Free Energy
Having discussed the Fokker-Planck equation that governs the evolution of the configurational distribution function, we have completed the system of equations for the unknown fields ρ s , v, θ s and ϕ. It remains to identify the fluxes, which are quantities that are unique for the given fluid (they provide a unique characterisation of the given material or given class of materials). The first step in this direction is the identification of energy storage mechanism, which is done by means of specific Helmholtz free energy.
We assume that the Helmholtz free energy ψ, [ψ] = J kg, is given by the formula where K is a constant to be specified later (so far its only purpose is to provide a normalisation factor in the argument of the logarithm function, since we must have a dimensionless quantity under the logarithm. We also recall the physical units of the configurational distribution function ϕ, [ϕ] = 1 m 6 , particle number density n p , n p = 1 m 3 , potentials U e , [U e ] = J, and U η , U η = J, and the Boltzmann constant k B , [k B ] = J K). The inclusion of the potentials U e and U η into the formula for the Helmholtz free energy is a standard one known from elasticity theory, while the term ϕ ln ϕ is a well known contribution due to configurational entropy of the dumbbells. The solvent contribution ψ s (θ s , ρ s ) to the specific Helmholtz free energy is left unspecified, since its particular form is not important for the ongoing discussion. We however give a particular formula for ψ s (θ s , ρ s ) in Section 6, where it plays a fundamental role.
Having introduced the formula for the Helmholtz free energy, we can find other thermodynamic quantities of interest via differentiation. Namely, the formula for the specific entropy η reads and the specific internal energy e is given by the formula

Constitutive Relations
Once we have specified the Helmholtz free energy we can proceed with the derivation of constitutive relations for the fluxes j e , j η , j ϕ,x , j ϕ,q , and the Cauchy stress tensor T. Our first objective is to identify the entropy flux j η and the entropy production ξ, such that the evolution equation for the specific entropy has the form The requirement on the nonnegativity of the entropy production then gives us a hint regarding the choice of constitutive relations for the fluxes j e , j η , j ϕ,x , j ϕ,q .
So far, we have identified an evolution equation for the internal energy (3c), and we know that the relation between the Helmholtz free energy and the internal energy and the entropy reads ψ(θ s , ρ s , ϕ) = e − θ s η. Differentiation of this relation yields Making use of (27) and the fact that η = − ∂ψ ∂θ s , we then see that (28) reduces to Now we need to manipulate (29) into the form (26). The lengthy but essential algebraic manipulations that allow us to achieve this objective are described in the following section, see Section 5.1. In particular, the desired final form of the entropy evolution equation is given in (51). The reader who is at the moment not interested in the necessary algebraic manipulations can go directly to Section 5.2, where we discuss the implications of the entropy evolution Equation (51) with respect to the choice of constitutive relations.

Evolution Equation for the Specific Entropy
Let us for a moment denote Using this notation the formula for the Helmholtz free energy reads and the evolution equation for the entropy (29) can be-with a slight abuse of notation-rewritten as Let us first consider the term ∂F ∂ϕ dϕ dt . Making use of the Fokker-Planck Equation (8), we see that We will first focus on the term B that contains the divergence operator with respect to the q variable. We see that The first integral vanishes by virtue of the Stokes theorem and the boundary condition (9) on ∂D. In the second term we first evaluate the gradient, and then we proceed, as follows Furthermore, the second term on the right-hand side can be manipulated, as follows (36) where we have again used the fact that the Stokes theorem and the boundary condition in configurational space (9) imply that the first term vanishes. As a matter of fact, we have used that ϕ ∂D = 0, which is a consequence of the particular choice of the flux j ϕ,q , see Section 3.4 for details. So far, we do not, however, have a formula for the flux j ϕ,q , and therefore we can not decide whether the term really vanishes. We will simply assume that the term vanishes, derive the formula for the flux j ϕ,q , and then we retrospectively check, that the derived formula for the flux j ϕ,q and boundary condition (9) are consistent with this assumption (this is indeed the case). We also note that the critical term can also be rewritten as (37) hence, if we wanted to keep it in the equations, we could easily do so. If we summarise the partial results, we therefore see that Next we focus on the term A in (33). We see that and, hence Now, we can finally return to the formula (32) for the material time derivative of entropy, and we find that which reduces to where we have used the evolution equation for the solvent density (3a). We divide (42) by the solvent temperature, and manipulate the energy flux term div x j e in the usual manner. We also split the term containing div x j ϕ,x into a part with a temperature-independent coefficient and a part with a temperature-dependent coefficient, and after some manipulations we get We introduce the notation and we split the Cauchy stress tensor T into the mean normal stress m and the traceless part T δ , where m = def 1 3 Tr T and T δ = def T − 1 3 (Tr T) I (this standard manipulation allows us to treat the volume-changing and volume-preserving deformations separately, see for example Málek and Průša [21] for details). We also do the same for the tensor P, and we rewrite (43) as Now, we are ready to manipulate the part with the temperature independent coefficient Substituting (47) back into (43) yields We rewrite the last equation as which helps us to identify the fluxes and entropy production terms. It remains to manipulate the term with the potential U e . We see that This manipulation implies that we can rewrite (49) in the final form as

Entropy Production and Constitutive Relations
Now, we are in a position to choose j η , j e , j ϕ,x , j ϕ,q , m and T δ in such a way that the right-hand side is nonnegative. We start with the flux j ϕ,x in physical space, the nonnegativity of the term − ∫ D k B ϕ ∇ x ϕ • j ϕ,x dq is granted if we make, for example, the simple choice where ζ denotes the hydrodynamic drag coefficient, (the presence of the hydrodynamic drag coefficient is motivated by the insight into the microscopic dynamics of the polymeric chains, see for example Barrett and Süli [31] for details). The choice (52) in turn implies that we have to fix the energy flux as We note that (53) is the same energy flux as in Öttinger and Grmela [18]. We also note that only the "energetic" part U e of the potential U enters the formula for the energy flux j e , which is an expected result.
Next, we fix the constitutive relation for the flux j ϕ,q in the configurational space; we set where the coefficient 2 ζ is again motivated by the microscopic dynamics of the polymeric chains. In particular, we want to keep the link between the Fokker-Planck equation and the associated stochastic process, see, for example, Barrett and Süli [31] for details. If we recall the relation between the force and the corresponding potential (12), we see that (54), in fact, reads where which, together with (52), gives upon substitution into (7) the Fokker-Planck equation in the form This is the specific Fokker-Planck equation that we have used in Section 3.
From the left-hand side of (51) we can infer the entropy flux j η . If we use the formula for j e , we see that the term under the div x operator reads where j ϕ,x is given by (52). We note that the term with the "energetic" part U e of the potential U has cancelled, and only the entropic part U η enters the entropic flux, which is an expected result. Finally, we fix the constitutive relation for the Cauchy stress tensor, namely for the mean normal stress m and the deviatoric part T δ . We need a choice that makes the term [T δ − P δ ] ∶ D δ + m − 1 3 Tr P + ρ 2 s ∂ψ s ∂ρ s + 2k B θ s n p div x v on the right-hand side of (51) nonnegative. If we want to model the solvent as a classical compressible Navier-Stokes fluid, then we set where ν andλ are the positive constants referred to as the shear viscosity and bulk viscosity. If we use (59), then the full formula for the Cauchy stress tensor reads whereλ = 3λ+2ν 3 . The first three terms on the right-hand side of (60) are the familiar ones. In particular, we recognise the formula for the thermodynamic pressure hence, we can rewrite (60) as The polymeric contribution T polymer to the Cauchy stress tensor T is given by the last two terms in (62). If we recall the definition of the tensor P, see (44), we see that This is the Kramers expression for the polymeric contribution to the Cauchy stress tensor. We note that the Kramers expression has not been assumed a priori: it is a consequence of the choice of energy storage and entropy production mechanisms. The same formula also follows from the GENERIC formalism, see Öttinger and Grmela [18], where it is also a consequence of the proposed modelling approach and it is not assumed a priori.

Temperature Evolution Equation
A cautious reader might notice that the energy flux j e and the entropy flux j η depend on the potentials U e and U η and on the normalisation constant K. However, if one shifts the value of a potential by a constant, then this operation should have no impact on the dynamics of the system. The same should hold for the constant K. Its particular choice should have no impact on the dynamics.
We show that this is indeed true. First, we observe that the shift of the value of the potentials and the choice of K have no impact on the balance of mass (3a) (evolution equation for the solvent density ρ s ) and on the balance of linear momentum (78b) (evolution equation for the velocity field v). Indeed, the Cauchy stress tensor, see (60), depends on the gradients of the potentials. Next we observe that the Fokker-Planck equation (evolution equation for ϕ) in the form (57) also contains only gradients of the potentials. It remains to check the evolution equation for the last quantity of interest, namely for the temperature.
In order to do so, we first need to explicitly formulate the evolution equation for the temperature. We start with the evolution equation for the entropy in the form (29), which is Because the entropy is given as the derivative of the Helmholtz free energy with respect to the temperature, η = − ∂ψ ∂θ s , we can use the chain rule on the left-hand side of (64), and with a slight abuse of notation we get Making use of (65) in (64) yields On the left-hand side, we have obtained the material time derivative of the temperature field θ s , and we can explicitly evaluate all terms on the right-hand side. We take into account the particular structure of the Helmholtz free energy (23), which is and we arrive at By virtue of the linearity of the polymeric part ψ p with respect to temperature, we get Next, we recall that the expression on the left-hand side of (69) is the classical definition of the specific heat at constant volume c V,s for the solvent only, we use the definition of the thermodynamic pressure p th,s for the solvent, see (61), the balance of mass (3a), and we get Now, we are in a position to use the formula for the energy flux j e , see (53), and the Fokker-Planck Equation (8) for the configurational distribution function ϕ, which yields Let us now focus on the last term. Integration by parts in the last term together with the boundary condition (9) and the constitutive relation (54) for j ϕ,q then gives us If we want to, we can also integrate by parts in the last term on the right-hand side, and we get where we have again exploited the fact that ϕ ∂D = 0. Now, we substitute back into (72), and we arrive at Finally, we use the formula for the Cauchy stress tensor (62), which is which, upon substitution into (75), yields This is the sought evolution equation for the temperature field. A brief inspection of (77) reveals that the evolution equation for the temperature is not affected by shifts of the values of the potentials U e and U η and the value of scaling constant K.

Summary
The system of governing equations for the solvent density ρ s (t, x), velocity v(t, x), solvent temperature θ s (t, x), and configurational distribution function ϕ(t, x, q) read, as follows: where the spring force F is given via the potentials U e and U η by the formula and the Cauchy stress tensor T is given by the formula The temperature evolution equation reads and the thermodynamic pressure p th,s and the specific heat at constant volume c V,s are given in terms of the derivative of the solvent specific Helmholtz free energy ψ s , as follows: We recall that the solvent density ρ s is in the given model identified with the density of the whole solvent-polymeric chains mixture. We note that if the spring potential only contains the entropic part, that is if U e = 0, then the governing equations simplify considerably. This system of governing equations is essentially the same as that reported by Öttinger and Grmela [18], albeit the model has been derived in a different manner. In fact our approach is similar to that used in the derivation of macroscopic viscoelastic rate-type models, see especially Rajagopal and Srinivasa [37], Málek et al. [38], Hron et al. [17], Málek et al. [39], Málek et al. [40], and Dostalík et al. [41]. For a coupling between microscopic and macroscopic models, see also Souček et al. [42].
The boundary condition in the configurational space D reads see discussion in Section 3.2. Because we are working with a model with centre-of-mass diffusion, we also need a boundary condition for the configurational distribution function in the physical space. If we consider a vessel Ω with an impermeable rigid wall, we set which, by virtue of (52), means that ∇ x ϕ • n x ∂Ω = 0. Regarding the remaining boundary conditions in physical space, we can opt for any standard boundary conditions. Because we have access to the velocity field v as well as to the Cauchy stress tensor T, we can choose from plethora of boundary conditions used for flows of polymeric fluids, see, for example, Hatzikiriakos [43] or Málek and Průša [21].

Stability
Following Coleman and Greenberg [44] and Coleman [19], see also Gurtin [45,46], Grmela and Öttinger [47] and Bulíček et al. [20] for further discussion, we can exploit the thermodynamic basis of the derived model in a rudimentary stability analysis of thermo-mechanical processes described by the corresponding governing equations. In particular, it is straightforward to investigate the finite amplitude (nonlinear) stability of the spatially homogeneous stationary states in a thermodynamically isolated container.
The spatially homogeneous equilibrium stationary state θ s,eq , ρ s,eq , ϕ eq and v eq is, in our case, given as θ s,eq =θ s , (81a) ρ s,eq =ρ s , (81b) The symbol n p is the spatially homogeneous polymer number density, which is n p = N Ω , where N is the total number of polymeric chains in the container and Ω is the volume of the container, the solvent densityρ s is the spatially homogeneous density, which isρ s = M Ω , where M is the total mass of the solvent in the container. The symbolθ s denotes a spatially homogeneous temperature field, and Mn p ,θ s stands for the equilibrium configurational distribution function (22) discussed in Section 3.5. Clearly, the quadruple (81) is a solution to the governing equations subject to boundary conditions that guarantee that the energy flux j e , j η , and j ϕ,x , see (53), (58), and (52), vanish on the container wall. We also have the standard boundary condition in the configurational space (9), which is where j ϕ,q is given by (54). If we denote then using the governing equations it is straightforward to see that, since all of the fluxes through the vessel wall vanish, the net total energy E tot is constant and the net total entropy S is a nondecreasing function in time. Moreover, the total mass of the solvent and the number of polymeric chains are also preserved in time.
Consequently, the following functional is a reasonable candidate for a Lyapunov type functional for the analysis of the stability of the stationary spatially homogeneous state with temperatureθ s , see, for example, Bulíček et al. [20] for details, where E tot denotes the net total energy at the spatially homogeneous stationary state (we use the nomenclature Lyapunov type functional, since we will only show the decay along trajectories and non-negativity of the functional everywhere except at the equilibrium. We will not investigate the link between the proposed functional and a suitable norm in the given state space, as this is beyond the scope of the current contribution). The rationale for the choice of such a functional is apparent from the following manipulation where ξ denotes the entropy production, which is a nonnegative quantity that vanishes at the spatially homogeneous stationary state (81). In the current case, we are forced however to modify the functional by explicitly adding zero terms which help us to articulate the constraints of mass conservation and the conservation of total number of polymeric chains. While the values of Lagrange multipliers can be found a priori, see, for example, Coleman [19] for a discussion in the case of a compressible fluid, we shall not proceed in that direction (in fact we have already tacitly identified one of the multipliers: the equilibrium temperatureθ s is the multiplier that enforces the conservation of net total energy). We shall find instead the multipliers on-the-fly. Indeed, from the pragmatic point of view, we will need the additional zero terms (after the integration) to show the nonnegativity of the integrands in V meq . Finally, we note that, in order to ensure that the functional vanishes at the equilibrium, we need to shift the functional by a constant value. We will denote the shifted functional by the same symbol, since this shift has no impact on the time evolution of the functional: it is just a matter of normalisation.
By virtue of (24) and (25), we can further split the internal energy and the entropy into their respective solvent part and polymeric part: e (θ s , ρ s , ϕ) = e s (θ s , ρ s ) + e p (θ s , ρ s , ϕ) , (89a) where e p (θ s , ρ s , ϕ) = def 1 Next, we split the function g into a solvent and a polymeric part as well g θ s , ρ s , ϕ θ s ,ρ s ,φ = g s θ s , ρ s θ s ,ρ s , + g p θ s , ρ s , ϕ θ s ,ρ s ,φ , where and, finally, we also split the whole functional, which is we write We shall first find an explicit formula for the polymeric contribution, see Section 6.2, and then we introduce an explicit equation of state for the solvent, and we find the formula for the solvent contribution, see Section 6.3.

Polymeric Part
First, we choose the characteristic temperature θ ref in the polymeric part of the Helmholtz free energy (23) asθ s , which is we set We are using the fact that the temperature fieldθ s at equilibrium is constant in space and time. Straightforward substitution of (90) into (92b) then yields Now, we can combine the terms with the potentials to deduce that Because the choice of the normalisation factor K has no impact on the dynamics of the system, we choose it to serve our objective. We fix the normalisation factor K as Using the notation introduced in (22), we see that (96) with the constant K chosen as in (97) reduces to on the closed system, namely, the constraint on the constant total number of polymeric chains. See the discussion in the introductory part of this section. We set V meq,p = Finally, using (106), we obtain the specific Helmholtz free energy ψ s defined by the thermodynamic relation ψ s (θ s , ρ s ) = e s (θ s , ρ s ) − θ s η s (θ s , ρ s ). The explicit formula for the specific Helmholtz free energy as a function of its natural variables reads and using this formula in (23) gives us a complete thermodynamic characterisation of the given fluid.
We are now in a position to construct the functional V meq,s , see (93). As in the previous section, we first choose the characteristic temperature of the solvent θ s,ref asθ s , which is we set We are using the fact that the temperature fieldθ s at equilibrium is constant in space and time. With this choice of the reference temperature field θ s,ref , we use the formulae for the specific internal energy e s , and the specific entropy η s , see (106), and we get Now, we need to manipulate the right-hand side of (109), and we need to show that it is nonnegative. This can be done, as follows. First, we exploit conservation of mass, which implies that ∫ Ω (ρ s −ρ s ) dv = 0. See also the discussion in the introductory part, in particular the discussion following (86). This observation allows for us to write ∫ Ω ρ s c V,sθs dv = ∫ Ωρ s c V,sθs dv. Second, we fix the reference pressure as and regroup the terms of (109) arriving at We note that the choice (110) leads to vanishing entropy at the equilibrium state. In this sense, the choice (110) is tantamount to the convenient fix of the entropy value at the equilibrium. Next we again exploit conservation of mass, which allows us to eliminate the last term on the right-hand side of (111), and yields Moreover, by adding the zero term to the right-hand side of (112), we arrive at Because the function f (x) = def x − 1 − ln x, is nonnegative for x ∈ (0, ∞) and vanishes if and only if x = 1, we obtain the desired nonnegativity of the first term on the right-hand side of (114).
Finally, we denote gρ s ,b (ρ s ) = def ρ s ln ρ ŝ ρ s where ρ s ∈ 0, 1 b . Clearly, gρ s ,b (ρ s ) = 0. Moreover, by calculating the first and the second derivatives of the function gρ s ,b one obtains that, irrespective of the values of b > 0 andρ s ∈ 0, 1 b , the function gρ s ,b is nonnegative and vanishes if and only if ρ s =ρ s . The second term on the right-hand side of (114) is thus nonnegative as well. We also note that the function gρ s ,b (ρ s ) blows-up as ρ s → 1 b −. Consequently, we have shown that the functional V meq,s is nonnegative and vanishes if and only if last term monitors the approach of the density to the equilibrium density. However, we recall that we have not discussed the relation between the functional and a metric structure on the state space; hence, the functional does not, at the moment, deserve to be referred to as the Lyapunov functional. See also Dostalík et al. [50,51] for a similar application in the case of macroscopic viscoelastic rate-type models, where the authors have actually found a relation between the proposed Lyapunov type functionals and a metric on the state space.
Inspecting (119) and (118), we also note that the time derivative of V meq is proportional to the net entropy production that vanishes at the stationary spatially homogeneous equilibrium state (81). These observations conclude our rudimentary stability analysis.
Finally, we note that, on the right-hand side of (119), we explicitly see all entropy production mechanisms that are active in the given fluid. Because the entropy production terms, in fact, determine the fluxes in the sense of the discussion in Section 5.2, we see that the identification of the entropy production and the identification of the Helmholtz free energy indeed provide a complete characterisation of the fluid of interest.

Conclusions
We have provided a simple derivation of a model for non-isothermal flows of dilute polymeric fluids. The analysis outlined above shows that if the microscopic dynamics of the polymeric chains is governed by the Fokker-Planck equation with a centre-of-mass diffusion term, then it is possible to find a corresponding evolution equation for the temperature field in such a way that the whole system of governing equations is compatible with the first and second law of thermodynamics.
However, we note that, although the derivation outlined above seems to be completely self-consistent, some of the steps require insight into the microscopic dynamics of the polymeric chains. In particular, the interpretation of the coefficient ζ as hydrodynamic drag, see especially Formula (55) for the flux j ϕ,q , would be impossible without the insight into the microscopic nature of the model (this is however true also for the derivation proposed by Öttinger and Grmela [18]). Without such an insight, one could easily misidentify the flux j ϕ,q and lose the connection between the Fokker-Planck type equation and the underlying stochastic process at the level of individual polymeric chains.
While the presence of the centre-of-mass diffusion term is motivated by physical considerations, see for example El-Kareh and Leal [52], Bhave et al. [53], Schieber [54] and Degond and Liu [32], it turns out that its presence can be gainfully exploited in the mathematical analysis of the corresponding governing equations in the purely mechanical setting as well, see especially Barrett and Süli [31,33,55,56,57] or Feireisl et al. [34]. (See also Bulíček et al. [58], Bulíček et al. [59] and Bathory et al. [60] in the context of macroscopic models with a stress diffusion term.) One might hope that once thermodynamically consistent models for the full thermo-mechanical coupling are derived, similar rigorous mathematical results might also be obtained for the corresponding full system of governing equations, that is including the governing equation for the temperature.
In particular, the knowledge of the thermodynamic background leads to natural a priori estimates, and to the design of natural functionals that might be used to monitor the stability of steady equilibrium or non-equilibrium states or to analyse weak-strong uniqueness. (See for example Feireisl and Novotný [61], Feireisl et al. [62] and Feireisl and Pražák [63] for the discussion in the case of a compressible Navier-Stokes-Fourier fluid. In the context of purely mechanical models of dilute polymeric fluids and stability of the corresponding flows, see especially Jourdain et al. [64].) Indeed, we formally-that is under the assumption that the governing equations posses a strong solution-identify such a functional, which might be of use in such analyses.
Author Contributions: All four authors were equally involved in the formulation of the problem. The calculations were carried out by M.D. and V.P. The discussion of the results and the presentation of the material were carried out by all four authors. All authors have read and agreed to the published version of the manuscript.