A Kinetic Perspective on k − e Turbulence Model and Corresponding Entropy Production

In this paper, we present an alternative derivation of the entropy production in turbulent flows, based on a formal analogy with the kinetic theory of rarefied gas. This analogy allows for proving that the celebrated k − e model for turbulent flows is nothing more than a set of coupled BGK (Bhatnagar–Gross–Krook)-like equations with a proper forcing. This opens a novel perspective on this model, which may help in sorting out the heuristic assumptions essential for its derivation, such as the balance between turbulent kinetic energy production and dissipation. The entropy production is an essential condition for the design and optimization of devices where turbulent flows are involved.


Introduction
Even though it is omnipresent around us, turbulence remains one of the last unsolved problems in classical physics.We can find turbulent flows in astrophysical plasmas, atmospheric winds, ocean currents, planetary boundary layers, rivers, fluid dynamics around all kinds of vehicles, internal combustion engines, heat exchangers and even blood flows.Turbulence indicates certain complex and unpredictable motions of a fluid with an entire hierarchy of structures called eddies.An overview of this topic is too broad to be reported here, and it can be found in [1].Interesting historical notes can be found in [2].In spite of the tremendous work done in the last century and the new possibilities offered by direct numerical simulation of turbulent flows, the theoretical progress in understanding turbulence is proceeding slowly.One of the reasons may be that turbulence is deeply rooted in fluid mechanics, while fluid dynamic equations are not derived from fundamental considerations and the mathematical theory of fluids is not yet fully developed [3].Up to now, the most important contribution to the theory of turbulence is still due to the Russian mathematician Andrei Kolmogorov.In 1941, a conceptual framework for turbulence (K41 theory) was defined by Kolmogorov [4][5][6].The K41 theory was originally derived for homogeneous and isotropic turbulence (i.e., statistically invariant under translations and rotations).After more than 70 years, it has been found out that the K41 theory may also explain turbulent flows over rough walls, by rationalizing the empirical scaling of Blasius and Strickler included in the Nikuradse's experiments [7].From a technological point of view, this finding has tremendous implications for heat transfer enhancement techniques [8][9][10].
Fundamental experiments in turbulence are extremely difficult because they are expensive and require long runs for collecting enough statistical data on the fluctuating quantities of interest.This is the reason why direct numerical simulations of turbulence are considered valuable alternatives.Von Neumann predicted in 1949 a revolution in the study of turbulence thanks to the advent of digital computers because of the possibility to simulate the fluid dynamic equations in 3D turbulent regimes [1].However, the first three-dimensional simulations of turbulence were carried out only many years later, when the advent of supercomputers and the development of (pseudo)-spectral numerical techniques based on fast Fourier transform algorithms unlocked adequate simulation capabilities: Steven Orszag and G.S. Patterson could achieve such a goal in 1972 [11].Clearly, numerical simulations offer the advantage of accessing the entire velocity field and other quantities of interest (e.g., vorticity, local energy dissipation, etc.) [1].The mesh size of these simulations (nowadays, up to tens of billions of grid-points) determines the spatial resolution, which depends on the maximum Reynolds number of the turbulent flow: the higher the Reynolds number, the bigger the mesh size.Due to recent progress in high performance computing, the feasible Reynolds numbers are approaching those recovered in experimental facilities.However, many times, history has shown that the main speedup in numerical computations comes from more efficient algorithms rather than better hardware [12,13].More than 20 years ago, the lattice Boltzmann method [14][15][16][17][18][19][20] was proposed as a powerful numerical method, which is prone to take advantage of the modern resources of high performance computing, especially for multi-physics and multi-scale simulations.The range of applications of the lattice Boltzmann method is very broad, including thermal radiation [21], thermal conduction [22], combustion [23][24][25][26], porous media [27,28], multi-component flows [29,30] and turbulence [31][32][33], to mention a few.In particular, it was argued that turbulence can be more efficiently modeled by the extended kinetic (Boltzmann) equation rather than continuum Navier-Stokes equations [31,32], by means of a proper discretization (lattice).Recently, a further improvement in computational performances was achieved by link-wise artificial compressibility method [34], which can be considered as a subsequent refinement of the lattice Boltzmann algorithm.Moreover, stability of kinetic schemes can be also improved using the pioneering idea of limiters [35].
In most of the engineering applications where turbulent flows are involved, direct numerical simulations are still unfeasible.Large eddy simulations and Reynolds averaged Navier-Stokes simulations are commonly used instead [36].Large eddy simulations using the lattice Boltzmann method have been carried out [37][38][39][40], as well as those using Reynolds averaged Navier-Stokes (RANS) equations [41,42].However, in the context of the lattice Boltzmann method, engineering turbulence models are still lacking, thus more research efforts are required [19].Among the closure models for the Reynolds averaged equations, the k − model is one of the most popular [36].The model takes its name from the turbulent kinetic energy k and the turbulent dissipation rate .The role of k − model in the theory of turbulence can be found in [43].The first two-equation model for predicting turbulent flows was proposed by Kolmogorov in 1942 [44], using fluctuation energy and frequency as its main variables.In 1968, Harlow and Nakayama proposed the k − model for turbulence [45], which was later refined [46].However, the true development of the model is credited to Launder, Hanjalic and Jones [47,48], whereas the collaboration with Spalding brought out the superiority of this model to a broad scientific audience [49][50][51].
In spite of its success from the practical point of view, the theoretical foundations of the k − model are still unclear.In fact, the model is made by two phenomenological equations and, in particular, the equation for the dissipation rate of turbulent kinetic energy is guessed by an analogy to the kinetic energy one [36].The idea to apply statistical mechanics to the study of fully developed turbulent flows is quite old, and it dates back to the pioneering work by Reynolds in 1883 [3].Indeed, the statistical approach is one of the main approaches to turbulence (along with the structural and deterministic ones).In particular, non-extensive statistical mechanics [52] and the recently proposed fluctuation theorems [53] appear promising in understanding turbulence.However, a fundamental problem still remains in investigating the relationship between the definition of entropy in case of non-equilibrium flows and the most popular turbulent models (e.g., k − model), which is exactly what we investigate here.
Recent trends in thermodynamics may help in achieving this goal.Starting from the theoretical framework introduced by Lars Onsager, Ilya Prigogine paved the way to the Thermodynamics of Irreversible Processes (TIP) by extending the fundamental Gibbs relation to non-equilibrium states [54,55].Such extension assumes that local equilibrium is achieved throughout the system, namely all subparts of the system are close to equilibrium condition and can be safely described by thermodynamic relationships and variables.Extended Irreversible Thermodynamics (EIT) [56,57] is another active field of investigation, where the local equilibrium assumption can be eventually overcome.The basic idea underlying EIT is that thermodynamic quantities are related to dissipative fluxes, which are considered as independent variables.Hence, the generalized Gibbs relation for non-equilibrium states can be obtained from the differential form of the generalized entropy, and it satisfies the following requirements: (1) in the limit of quasi-static processes, it can be reduced to the equilibrium Gibbs relation; (2) it is Galilean invariant; (3) results are supported by experimental evidence [58].This procedure is not even limited to bulk flows, as it can be successfully adopted for studying non-equilibrium interfaces [59][60][61].
Surprisingly enough, only recently Extended Irreversible Thermodynamics has paid attention to turbulent flows [62][63][64][65].There are two main approaches in literature.In the first approach, Jansen [66] and Hauke [67] have extended entropy-based stability analysis to turbulent flows.The key idea consists in applying the Reynolds time averaging to the entropy transport equation obtained by the Thermodynamics of Irreversible Processes (TIP), i.e., by assuming that local equilibrium assumption holds.This approach has been somehow re-proposed recently and even applied to commercial codes [68][69][70][71].However, the entropy transport equation is highly non-linear, and several ways to perform the Reynolds time averaging could be adopted.For example, multiplying this equation by the temperature before applying the time-averaging leads to simplified calculations, but a closure problem appears (for fluctuating temperature and entropy production correlations) and different models have been proposed [72].In fact, the mean entropy production cannot be easily expressed by other mean flow variables alone [72].We believe that these difficulties arise from understanding the local equilibrium assumption.Hence, we aim to clarify the derivation of entropy production in turbulent flows.The entropy production is helpful in the design and optimization of devices dealing with turbulent flows (e.g., by the entropy generation minimization approach [73,74]).
In this work, materials and methods representing the starting point of our analysis are briefly summarized in Section 2; then, the main results for turbulent flows are reported in Section 3, where the k − model is recast by a kinetic approach.Finally, Section 4 presents some considerations derived from the reported results; whereas conclusions are drawn in Section 5.

Materials and Methods
In the asymptotic limit of low Mach number flows, the incompressible limit of the Navier-Stokes-Fourier system of equations for Newtonian fluids with constant transport coefficients [75,76] can be written as where ρ 0 is the average fluid density (assumed constant and different from the actual fluid density ρ), u is the fluid velocity, p is the pressure, ν is the kinematic viscosity, h is the enthalpy, λ is the thermal conductivity, and ∇ S u = 1/2 ∇u + ∇ T u is the strain rate tensor.Alternative formulations of Equation ( 1) may be expressed in terms of the internal energy or total energy, instead of enthalpy.
It is well known that, in regimes where momentum convection prevails on momentum diffusion (u • ∇u ν∇ 2 u) or, equivalently, Reynolds number is larger than a flow-dependent threshold, the flow field is characterized by rapid variations of pressure and velocity in space and time, i.e., by chaotic coherent structures called turbulent eddies.This discussion can be more rigorous by adopting characteristic quantities, which can be defined in different ways according to what they refer to (fields or operator).On the one hand, the characteristic quantity u c gives the magnitude of the velocity field, being u c = max(u) in the considered domain.On the other hand, operators may involve more difficult definitions.For instance, the characteristic length l c is defined by where ϕ is a generic field property of the fluid flow.Hence, ).This clarifies the reason why the inequality u • ∇u ν∇ 2 u implies large Reynolds numbers: Re = u c l c /ν 1.As originally suggested by Reynolds in 1895, let us consider the time-averaged equations of motion for fluid flows, i.e., the so-called Reynolds-averaged Navier-Stokes equations (RANS) [36].First, let us introduce the filter operator • for the generic quantity where ϕ is assumed independent on the initial condition t (in the context of chaotic dynamical system, this means that the system can only have one strange attractor).
From a practical point of view, the limit in Equation ( 5) is truncated to some finite time, which is assumed to be much larger than the characteristic time of turbulent fluctuations.The previous definition allows one to introduce the average velocity field ū = u and the corresponding velocity fluctuation u = u − ū or, equivalently, to decompose the actual velocity field, namely u = ū + u .Introducing the previous decomposition in the divergence-free condition given by Equation ( 1) and taking the time average • shown in Equation ( 5) yield and consequently ∇ • u = 0, meaning that both the average and fluctuating fields are divergence-free.Proceeding in a similar way for Equation (2) yields where p = p .The term u • ∇u represents the novel physical content of the RANS momentum equation.Of course, in order to derive a set of closed RANS equations, the latter term must be modeled by quantities depending on the averaged quantities (or their gradients) only.First of all, this term must be made symmetric.Taking into account that ∇ • u = 0, the following condition holds Next, let us consider the Boussinesq's eddy viscosity assumption, namely where ν t is the eddy viscosity (additional kinetic viscosity due to turbulent eddies beyond molecular one) and ∇ S ū is the strain rate tensor defined as ∇ S ū = (∇ ū + ∇ ūT )/2.It is easy to prove that ∇ ū = ∇ S ū + ∇ W ū and, more importantly, ∇ S ū : ∇ W ū = 0, where ∇ W ū = (∇ ū − ∇ ūT )/2 is the vorticity tensor.Substituting the Boussinesq's eddy viscosity assumption into Equation (7), and taking into account that 2 ∇ • ∇ S ū = ∇ 2 ū (because of Equation ( 6)), yield It is clear that working with RANS is a successful approach, as far as an accurate model is defined to compute the corrective factor (ν + ν t ) in Equation (10).Clearly, the key idea is to investigate the part of kinetic energy due to turbulent fluctuations.Let us substitute ∇ 2 u with 2 ∇ • ∇ S u in Equation ( 2) and let us multiply the result by u, namely where e k = u 2 /2 is the total kinetic energy.After some simple algebra, the previous equation can be rewritten as where p k = p/ρ 0 is the kinetic pressure and k = 2ν(∇ S u) 2 is the (positively defined) dissipation function of the total kinetic energy (the square power (∇ S u) 2 simply means ∇ S u : ∇ S u).Similarly to what has been done for the velocity field, let us now decompose both the total kinetic energy e k and its dissipation function k , namely e k = e k + e k and k = k + k , respectively.In this case, the main difference is that the time average of the previous quantities, e.g., e k , does not coincide with the same definition applied to averaged quantities, e.g., ēk .In particular, the following relations hold where ēk = ū2 /2, k = (u ) 2 /2 is the turbulent kinetic energy, ¯ k = 2ν(∇ S ū) 2 and = 2ν (∇ S u ) 2 is the turbulent dissipation function.For example, the previous definitions imply that e k = ū • u + u • u /2 − k and consequently e k = 0, even though e k = ēk .The turbulent kinetic energy k is the excess of kinetic energy due to turbulent fluctuations beyond the one caused by the average flow.Similarly, the turbulent dissipation is the excess of dissipation due to turbulent fluctuations.These are the two main quantities defining this approach to turbulence modeling, and this is the reason why this method is called k − model [36].The second key idea is that the eddy viscosity ν t depends on k and only.The physical dimensions of are those of kinetic energy divided by time, thus it can be interpreted as the dissipation rate of turbulent kinetic energy.Based on dimensional analysis, the following phenomenological relation is assumed for the eddy viscosity where C t is a tunable constant of the model (see next).The next step consists in deriving a closed set of equations for k and .Let us substitute the previous decompositions in Equation (12) and then perform the time averaging.Some difficulties may arise from cubic (third-order) terms with respect to velocity, i.e., e k u , but the following condition allows one to sort them out: where the last term can be expressed by the Boussinesq's eddy viscosity assumption Equation (9).Exploiting the condition in Equation ( 16), it follows: On the other hand, multiplying Equation ( 7) by the average velocity ū yields Subtracting Equation (18) from Equation ( 17) and taking into account that k = e k − ēk (because of Equation ( 13)) yield Taking into account the following equivalence Equation ( 19) can be rewritten as which is the standard form for the turbulent kinetic energy equation (TKE).Some simplifications are usually applied to the flux of turbulent kinetic energy.First of all, the last term of such a flux is usually negligible at high Reynolds number flows, because it is proportional to νν t (even though ν t can be up to roughly ten times ν, still this term is proportional to ν 2 ).Secondly, the (leading) term (u • u /2) u clearly shows the turbulent transport of the quantity u • u /2, which is the argument of the turbulent kinetic energy k.Hence, the key idea is to generalize the gradient-diffusion approximation [36], namely where σ k is a tunable constant of the model (see next).Consequently, the equation of the turbulent kinetic energy in the k − model becomes In the latter equation, the sign of the right hand side is no more uniquely prescribed.Here, the quantity (ν t /ν) ¯ k acts as turbulent energy production, moving kinetic energy from the mean flow to the turbulent fluctuations; while acts as turbulent energy dissipation, moving energy in the opposite direction.The same theoretical framework is used to derive the equation for turbulent dissipation.However, the equation of turbulent dissipation in the k − model is completely heuristic, being substantially derived by analogy with the previous equation [36].The source/sink of turbulent dissipation is derived by dividing the right hand side of Equation ( 23) by a proper characteristic time (∼ k/ ) and introducing ad hoc some tunable constants, namely where σ , C 1 and C 2 are tunable constants of this model.The standard k − model is defined by Equations ( 6) and (10) for the average fluid flow, by Equation ( 15) for the eddy diffusivity, by Equations ( 23) and ( 24) for the turbulent kinetic energy and turbulent dissipation and, finally, by the following set of constants, C t = 0.09, C 1 = 1.44,C 2 = 1.92, σ k = 1.0, σ = 1.3 [36].

Results
Let us start considering the entropy production of (incompressible) laminar flows.Non-equilibrium thermodynamics of laminar flows can be obtained by an entropy production equation, which can be derived by following the typical guidelines of TIPs [54,55].Let us assume that all subparts of the system are close to equilibrium conditions, so that they can be described by classical thermodynamics.In equilibrium conditions, the fundamental Gibbs relation allows for relating entropy with other thermodynamic potentials, namely Tds = dh − dp/ρ 0 , being T the temperature and s the entropy.Here, we assume that this relation holds also in non-equilibrium conditions, where the Gibbs relation becomes the entropy definition.However, this definition must be Galilean invariant and thus Lagrangian time derivatives must be considered, namely Substituting Equation (3) into Equation ( 25) yields The previous equation can be rewritten as where α = λ/(ρ 0 c p ) is the thermal diffusivity.By further elaborating on the first term at the right hand side of the previous expression yields or equivalently where Note that the prime notation is used to indicate the entropy production per unit of mass, i.e., σ , instead of the more common production per unit of volume, i.e., σ = ρ 0 σ .In the previous equations, each specific transport phenomenon is a source of entropy production, namely: σ α and σ ν represent the entropy produced by heat transfer (ruled by temperature gradient) and by fluid flow (ruled by strain rate), respectively.Clearly, each entropy source is in agreement with the second law of thermodynamics, namely σ α ≥ 0 and σ ν ≥ 0. The same considerations can be made for the global entropy production, namely σ α + σ ν ≥ 0. This confirms that the hypothesis reported in Equation ( 25) is valid for non-equilibrium states, at least as far as laminar flows are concerned.
In case of turbulent flows, the situation is more complex.First, the entropy production due to velocity gradients σ ν (Equation ( 31)) bears a strong resemblance with the dissipation function k appearing in the kinetic energy Equation ( 12), namely The above observation paves the way to a meaningful formal analogy.In fact, kinetic models are often used to explain microscopically the entropy production; hence, the previous observation suggests that the same kinetic models may be useful to explain also the turbulent dissipation function k and, consequently, the turbulent kinetic energy k.The simplest kinetic model is the celebrated Bhatnagar-Gross-Krook (BGK) model [77].Here, we propose interpreting the k − turbulence model as a set of coupled BGK-like equations.Let us reformulate Equation (23), which is the only one rigorously derived in k − model, as follows: where τ k = k/ is the characteristic relaxation time and k eq is the local equilibrium for the turbulent kinetic energy, namely Similarly, Equation ( 24) can be recast as where τ = τ k /C 2 is the second characteristic relaxation time and eq is the local equilibrium for the turbulent dissipation function, namely The last term in Equation ( 35) is an empirical forcing.It is worth the effort to elaborate further on the analogy with the kinetic equations for rarefied gas.In the Boltzmann equation, the external forcing is described by the term g • ∇ v f , where g is the acceleration vector due to the external force field, v is the single-particle velocity vector and f is the single-particle distribution function [78].Assuming f ≈ f eq , where f eq is the Gaussian distribution function, yields g • ∇ v f ≈ −g • (v − u)/e f eq , where e is the specific internal energy.In the last expression, g • (v − u)/e is the characteristic frequency of the external force field, which remarkably corresponds to 1/τ (1 − C 1 /C 2 ) in the last term of Equation (35).Hence, the celebrated k − model for turbulent flows is nothing more than a set of two coupled BGK-like Equations ( 34) and (36) with an empirical forcing.Realizing this formal analogy represents the most important result of the present paper, with important consequences discussed in the following.
Before proceeding further, it is worth highlighting an important simplification.In Equation ( 28), some terms proportional to the product between (1/T) (generalized intensive quantity) and ∇T (generalized force) appear.In the following, Reynolds decomposition will be considered only if the temperature is the argument of a spatial gradient (owing to the small size of the turbulent eddies).On the other hand, the average value of temperature will be considered for the intensive quantity 1/ T, because turbulent fluctuations around the room temperature are negligible (this is not the case for u , because ū can be eventually zero).Introducing the usual Reynolds decompositions for T (in generalized force only), s, h and u, and performing the time averaging • yield Let us consider the advection term, namely u • ∇s = ∇ • s u .The argument of the divergence operator can be expressed by the gradient-diffusion approximation [36], namely where β and γ are tunable transport coefficients and, in particular, Substituting Equation (38) into Equation (37) yields We now have to find out some approximated expressions for the main quantities, namely k and .Let us consider again the BGK-like Equation (33), which is more rigorous than Equation (35) and is also easier to be analyzed because there is no forcing term.The main fluid flow characteristic time is defined as τ c = l c /u c (see above).Clearly, the characteristic relaxation time of turbulent kinetic energy τ k is much shorter than the fluid flow characteristic time, namely τ k /τ c 1, because the turbulent structures are much smaller.Hence, it is possible to find out an approximated solution of the set of BGK-like equations defining the k − model in the asymptotic limit τ k /τ c 1.For our purposes, it is enough to consider the following formal expansion, without expanding the differential operators [78], namely Substituting the previous ansatz into Equation (33) and collecting terms with the same order of magnitude with regards to τ k /τ c [78], we obtain k (0) = k eq at the leading order.This means that k = k eq + O(τ k /τ c ) and k ≈ k eq in the asymptotic limit τ k /τ c 1 or, equivalently, ≈ eq .Substituting ≈ eq into Equation (40) yields The derivation proposed here is based on the formal analogy of Equation ( 33) with kinetic equations, in particular with BGK-like equations.However, this is consistent with the canonical derivation based on ≈ (ν t /ν) ¯ k , which means that there is substantially a balance between turbulent kinetic energy production and dissipation.The condition = (ν t /ν) ¯ k also implies which allows computing (∇ S u ) 2 by means of (∇ S ū) 2 , where the latter term is usually the only one available in practical calculations.The same idea can be applied to simplify the entropy production due to temperature gradient as well.It is easy to verify that the balance between turbulence production and dissipation implies where α t = ν t /Pr t and Pr t is the turbulent Prandtl number.Substituting Equation (44) into Equation (42) yields which is valid for both laminar and turbulent flows, and it represents the generalization of the second law of thermodynamics for turbulent flows.

Discussion
Equation ( 45) confirms the applicability of the Thermodynamics of Irreversible Processes to (incompressible) laminar and turbulent flows.Some comments are reported in the following.
1.In this paper, we present an alternative derivation of Equation ( 45) based on a kinetic approach.
This approach provides a novel perspective on the k − turbulence model, which remains one of the most successful models for many engineering applications, even though it is still affected by empirical assumptions for the turbulent dissipation function.In particular, our approach clarifies that this model is nothing more than a set of coupled BGK-like equations with a proper forcing, see the last term in Equation (35).Note that the formal expansion proposed for k in Equation ( 41) may not be suitable for , because the forcing term would change the local equilibrium, which must be unique in kinetic theory.Hence, some further investigations are required to find out the most suitable expansion for analyzing the high-order asymptotics.This is not surprising, because different asymptotic approaches (Chapman-Enskog, Hilbert, Grad, etc.) have been long debated in kinetic theory of rarefied gas [78].2. Equation (45) proves that four terms are the main sources of entropy production rates in turbulent flows: (1) direct dissipation; (2) indirect (turbulent) dissipation; (3) heat conduction driven by average temperature gradients and (4) heat conduction driven by fluctuating temperature gradients [68][69][70].These production rates (and their sum) are positively defined, consistently with the second law of thermodynamics.3.In Equation ( 45), the entropy production can be expressed as being η k the phenomenological coefficients of irreversible phenomena (η 1 = (α + α t ) c p / T2 and η 2 = (ν + ν t ) 2/ T) and X k the generalized thermodynamic forces (X 1 = ∇ T and X 2 = ∇ S ū).We selected k subscript as equal to the dimensionality of the corresponding thermodynamic force (vector k = 1, tensor k = 2).Consequently, X 2 k = X k * X k , where the generalized product * means scalar product • for k = 1 and saturation product : for k = 2.

Conclusions
In this paper, we present an alternative derivation of the entropy production in turbulent flows, based on a formal analogy with the kinetic theory of rarefied gas.This analogy allows for proving that the celebrated k − model for turbulent flows is nothing more than a set of coupled BGK-like equations with a proper forcing.This opens a novel perspective on this model, which may help in sorting out the heuristic assumptions essential for its derivation: for example, the balance between turbulent kinetic energy production and dissipation.Moreover, this has also important implications for applications, because the second law of thermodynamics for turbulent flows given by Equation ( 45) is very popular in the design and optimization of devices dealing with hypersonic flows (e.g., by the entropy generation minimization approach [73,74]).The new perspectives introduced by this work may find application in a broad variety of fields, spanning from aerospace (e.g., aeroelastic study of flexible flapping wings [79]) to materials (e.g., nature-inspired structures for fluid drag reduction [80]) research, from biomedical (e.g., turbulent blood flow [81]) to heat transfer (e.g., turbulent heat transfer in heat exchangers [82]) applications, from automotive (e.g., air conditioning components [83]) to atmosphere [84] modeling.