Irreversible Thermodynamics and Modelling of Random Media

1. Introduction

Engineering and Nature very often are concerned with media that do contain some randomness. One particularly well known example is given by turbulent flows, where random fluctuations are the results of the growth of unstable motions from small perturbations in the initial or boundary conditions, when the velocity gradients or/and temperature gradients are too large in some place within the flow. These fluctuating motions result in effective additional friction, effective additional diffusion of species or heat, and additional energy dissipation (i.e. transfer from organized kinetic energy to turbulent kinetic energy, and finally to internal energy). Another example of random flows is seen in two-phase flows, either with particles dispersed into a liquid or gaseous continuous phase, or built with a bubbly liquid. Here, a first reason of randomness is given by the fact that the locations and velocities of the particles (or bubbles) at initial time and in entrance sections cannot be known or controlled, and the induced fluctuations are not rapidly damped with time. Generally, this randomness becomes rapidly uncorrelated with these initial conditions due to the complex interactions between the flow of the continuous phase and the moving inclusions, interactions that modify also the number and sizes of these inclusions. Consequently, the flow conditions in both phases are or become turbulent, with interconnected fluctuations of velocities. Again, the effective result of the fluctuations is additional friction, additional diffusion or dispersion, and additional energy dissipation. Two-phase or multiphase flows with a large amount of solid particles are used in the industrial devices called as “Fluidized Beds”, and are encountered naturally in “Granular Flows”, and in this case there are many lasting contacts between the particles. The so called “Granular Media”, even without actual flow but just with some deformations or slow motions, are also subjected to randomness, due to the preparation of the medium but also to the differences between the grains shapes, and even the global properties of these media are very difficult to predict.

The properties of these kinds of media and flows have been studied since more than hundred years, and very useful prediction methods have been proposed so far, with methodologies becoming more and more clearly similar. Although all these media are showing clearly irreversible and dissipative processes, the typical reasoning of irreversible thermodynamics has never been used, at least explicitly (it might be present, but not consciously, in the brain of the authors). Our purpose here is to show that the Irreversible Thermodynamics and the “Second Principle” can justify many features of the models that have been proposed until now for these three types of applications. It has to be emphasized that Irreversible Thermodynamics deals with the modelling of media in evolution, whatever can be this model. The usual model of classical continuous medium, where the randomness lies only at the infinitely small scale of moving molecules, is not the only one that can be addressed. The so called “Extended Irreversible Thermodynamics” (Jou et al., 2001) allows us indeed to build pertinent entropy functions adapted to each type of model, in particular depending of the basic variables of the model used, as we will see.

The first part of this paper deals with classical turbulent flows. The general problem of building a model for these flows, i.e. an approximate mathematical representation that could give predictions close to experiments, has in practice given rise to several solutions, which are valid in different domains, i.e. for different kinds of experimental situations. More precisely, the domains of efficiency of the different models appear to be larger and larger for models using an increasing number of variables. The basic variables that allow the building of these models are “averaged variables” (now precisely defined as statistical mean values), in a finite number, and then they do clearly constitute a truncated representation of the medium, but it is expected that these variables are sufficient for the purpose of the approximated knowledge needed. Some of these variables do satisfy balance equations, which are obtained averaging the primitive balance equations of the mechanics of continuous media, but these equations do need “closure assumptions” before to be useful in constituting a “Turbulence Model”. Once the basic variables are chosen, the “Extended Irreversible Thermodynamics” allow us indeed to define an entropy function for this representation of the medium, for each of the models proposed. We will see that this entropy function is not the statistical mean value of the classical entropy for the fluid used in this turbulent flow, but has to be built taking into account the averaged form of the state equations of the fluid. Then, the Second Principle may lead to conditions concerning the ingredients of the turbulence model, conditions that are different for the different models. We will show that here simply considering two very popular models, and we will see that the usual practice does satisfy these conditions.

The second part of the chapter considers the modelling of two-phase flows. This modelling can be attacked with a similar approach and similar tools as turbulent flows: a set of averaged variables are defined; balance equations are written from the primitive equations of a piecewise continuous medium and closure assumptions are needed, to be adapted to each kind of two-phase flows. Here the closure assumptions are needed for both additional dispersion fluxes and exchange terms between the phases. These two-phase flows are again associated with irreversible effects, and again the framework of extended irreversible thermodynamics helps the modelling. We will recover again assumptions that are used in practice and supported by experiments.

The third part is devoted to granular media, with the same point of view. The approach generalizes the case of two-phase flows, because the granular medium is considered as a multiphase medium where each grain is considered as one phase, and the fluid in between as the last phase (Borghi § Bonelli, 2007). Again, the variables of the model are defined as statistically averaged grains phase and fluid phase variables, and their balance equations can be found from the primitive equations of mechanics, to be completed with closure assumptions. The approach can handle the situations where the medium is “quasi-static” as well as the ones where the medium is flowing with large velocities, giving rise to the limiting case of “Granular Gases”. In both cases, Extended Irreversible Thermodynamics can be used for building entropies and giving closure assumptions, without linear laws in this case particularly. The classical models of “Granular Gases” are recovered in the limit of large velocities, with the difference that the grain size is not the single length scale to be considered. For quasi-static situations, it is necessary to include in the model the six components of the averaged “Contact Cauchy Stress Tensor” of the grains phase. Intermediate situations can be handled as well.

2. Irreversible thermodynamics for turbulent flows models

Turbulent flows at present time are studied with models of two different types, namely Reynolds averaged Navier-Stokes (RANS) models, or Large Eddies Simulations (LES) models. We consider here for simplicity the framework of RANS models, although the extension to LES models may be not difficult. Then, the object of the study is not one given flowing medium, but an infinite number of such media, flowing submitted to initial and boundary conditions randomly perturbed. Indeed, what we call “one turbulent flow” is this ensemble of similar individual flows, because only statistical mean values calculated from this ensemble (or, at best, probability density functions) can be possibly predicted by some model.

The “mean medium” that is the subject of modelling is then described by “mean velocities”, “mean densities”, “mean enthalpy per unit of mass”, etc.. For N realizations, one defines:

$$\overline{\rho }=\lim (N\to \infty )(\frac{1}{N}{\displaystyle \sum _{n=1}^N{\rho }_{}^n})$$E1

at each given location within the flow field, and similarly for pressure. But in the general case where the density of the fluid can vary noticeably, it has been found more interesting to consider for other variables the mean value weighted by density, and one defines:

$${\tilde{v}}^{\alpha }=\lim (N\to \infty )({\displaystyle \sum _{n=1}^N{\rho }_{}^n}{v}_{}^{\alpha n}/{\displaystyle \sum _{n=1}^N{\rho }_{}^n})$$E2

,

i.e. the mean velocity is nothing but the mean momentum divided by the mean density. Similarly for the internal energy per unit mass (or the enthalpy per unit mass), the mean values are weighted with density.

Within a turbulent flow, these quantities are defined for each location and for each time, and vary with the location and eventually with time, but these variations are much smoother than the local and instantaneous quantities measured in one single realization of the flow, and by definition they are not sensitive to perturbations of initial and boundary conditions, that are inevitably present.

Is it possible to characterize the turbulent medium by $$\overline{\rho },\tilde{e},{\tilde{v}}^{\alpha }$$only, when a non turbulent flow of the same fluid is correctly characterized by$$\rho ,e,v^{\alpha }$$? Of course, the random fluctuations do have some influence. However, in some situations, it may be possible to take into account this influence by some algebraic formulas in terms of the mean variables and their derivatives. Of course, these formulas are only approximations, with a limited domain of validity assessed by experiments.

The first invented turbulence model, called “Prandtl mixing length model”, is of this type, and has been very useful in many practical cases. But a larger domain of validity can be covered taking into account additional fields variables pertinently chosen, described with their own partial differential evolution equations. There are different turbulence models of this improved type, the simplest one introduces only as additional field variable the kinetic energy of the random turbulent fluctuations, say$$\tilde{k}$$, the “length scale” characterizing the two-point correlation of the fluctuations, $$l_t$$, being given with algebraic formulas. More detailed models use $$\tilde{k}$$and$$\epsilon$$, the dissipation rate of$$\tilde{k}$$, described by partial differential evolution equations, or $$\tilde{k}$$and$$\omega$$, a characteristic frequency of turbulence, or the Reynolds stress tensor components themselves and$$\epsilon$$(Schiestel, 2006).

2.1. Basic equations

The physical basic equations corresponding to balance of mass, momentum and total energy give birth in any case to balance equations for the mean variables$$\overline{\rho },\tilde{e},{\tilde{v}}^{\alpha }$$. We write them with a tensor notation and the Einstein convention of repeated indices. We assume an orthogonal reference frame and then upper or lower indices are equivalent:

$$\frac{\partial }{\partial t}\left(\overline{\rho }\right)+\frac{\partial }{\partial x_{\alpha }^{}}(\overline{\rho }{\tilde{v}}_{}^{\alpha })=0$$E3

$$\frac{\partial }{\partial t}\left(\overline{\rho }{\tilde{v}}^{\beta }\right)+\frac{\partial }{\partial x_{\alpha }^{}}(\overline{\rho }{\tilde{v}}_{}^{\alpha }{\tilde{v}}^{\beta })=\frac{\partial }{\partial x_{\alpha }}(-\overline{p}{\delta }^{\alpha \beta }+{\overline{\tau }}^{\alpha \beta }-\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }})+\overline{\rho }{g}_{\beta }$$E4

$$\frac{\partial }{\partial t}\left(\overline{\rho }{\tilde{e}}_t\right)+\frac{\partial }{\partial x_{\alpha }^{}}(\overline{\rho }{\tilde{v}}_{}^{\alpha }{\tilde{e}}_t)=\frac{\partial }{\partial x_{\alpha }}(-{\overline{j}}_{\alpha }^Q-\overline{p{v}^{\alpha }}+\overline{{\tau }^{\alpha \beta }{v}^{\beta }}-\overline{\rho v{\text{'}}^{\alpha }{e}_t})+\overline{\rho }{\tilde{v}}^{\beta }{g}_{\beta }$$E5

The total energy here comprises the internal energy, the kinetic energy of the mean motion, and the kinetic energy of turbulent fluctuations, so$$\overline{\rho }{\tilde{e}}_t=\overline{\rho }\tilde{e}+\frac{1}{2}\overline{\rho }{\tilde{v}}^{\alpha }{\tilde{v}}^{\alpha }+\frac{1}{2}\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\alpha }}$$. We have defined$$v{\text{'}}^{\alpha }=v^{\alpha }-{\tilde{v}}^{\alpha }$$, and so$$\overline{\rho v{\text{'}}^{\alpha }}=0$$. The well known « Reynolds stress tensor » is $$R_{\alpha \beta }$$such that:$$\overline{\rho }{R}_{\alpha \beta }=\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }}$$. Similarly$$e\text{'}=e-\tilde{e}$$, but$$p\text{'}=p-\overline{p}$$. The “turbulent kinetic energy” introduced above is nothing but:$$p\text{'}=p-\overline{p}$$.

From (2) and (1), one obtains the equation for the kinetic energy of the mean motion:

$$\begin{array}{l}\frac{\partial }{\partial t}\left(\overline{\rho }{\tilde{v}}^{\beta }{\tilde{v}}^{\beta }/2\right)+\frac{\partial }{\partial x_{\alpha }^{}}(\overline{\rho }{\tilde{v}}_{}^{\alpha }{\tilde{v}}^{\beta }{\tilde{v}}^{\beta }/2)=\frac{\partial }{\partial x_{\alpha }}(-\overline{p}{\tilde{v}}^{\alpha }+{\tilde{v}}^{\beta }{\overline{\tau }}^{\alpha \beta }-\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }}{\tilde{v}}^{\beta })\\ +(\overline{p}{\delta }^{\alpha \beta }-{\overline{\tau }}^{\alpha \beta }+\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }})\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}+\overline{\rho }{\tilde{v}}^{\beta }{g}_{\beta }\end{array}$$E6

Then, from (3) and (4) we have:

$$\begin{array}{l}\frac{\partial }{\partial t}\left(\overline{\rho }(\tilde{e}+\tilde{k})\right)+\frac{\partial }{\partial x_{\alpha }^{}}(\overline{\rho }{\tilde{v}}_{}^{\alpha }(\tilde{e}+\tilde{k}))=\frac{\partial }{\partial x_{\alpha }}(-{\overline{j}}_{\alpha }^Q-\overline{pv{\text{'}}^{\alpha }}+\overline{{\tau }^{\alpha \beta }v{\text{'}}^{\beta }}-\overline{\rho v{\text{'}}^{\alpha }e}-\frac{1}{2}\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }v{\text{'}}^{\beta }})\\ -\overline{p}\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\alpha }}+({\overline{\tau }}^{\alpha \beta }-\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }})\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}\end{array}$$E7

It is possible also to derive directly from the non-averaged momentum equation a balance equation for$$\overline{\rho }\tilde{k}=\frac{1}{2}\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\alpha }}$$. That is done by scalar multiplying the instantaneous (non-averaged) momentum equation by the velocity in order to obtain an equation for the total kinetic energy, and then by subtracting the equation of the kinetic energy of the mean motion (4). One gets:

$$\frac{\partial }{\partial t}(\overline{\rho }\tilde{k})+\frac{\partial }{\partial x_{\alpha }^{}}(\overline{\rho }{\tilde{v}}_{}^{\alpha }\tilde{k})=\frac{\partial }{\partial x_{\alpha }}(\overline{v{\text{'}}^{\beta }{\tau }^{\alpha \beta }}-\overline{v{\text{'}}^{\alpha }p}-\frac{1}{2}\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }v{\text{'}}^{\beta }})-\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }}\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}+\overline{p\frac{\partial v{\text{'}}^{\alpha }}{\partial x_{\alpha }}}-\overline{{\tau }^{\alpha \beta }\frac{\partial v{\text{'}}^{\beta }}{\partial x_{\alpha }}}$$E8

From (5) and (6), finally the equation for the internal energy is written:

$$\frac{\partial }{\partial t}\left(\overline{\rho }\tilde{e}\right)+\frac{\partial }{\partial x_{\alpha }^{}}(\overline{\rho }{\tilde{v}}_{}^{\alpha }\tilde{e})=\frac{\partial }{\partial x_{\alpha }}(-{\overline{j}}_{\alpha }^Q-\overline{\rho v{\text{'}}^{\alpha }e})-\overline{p}\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\alpha }}-\overline{p\frac{\partial v{\text{'}}^{\alpha }}{\partial x_{\alpha }}}+\overline{{\tau }^{\alpha \beta }\frac{\partial v{\text{'}}^{\beta }}{\partial x_{\alpha }}}+{\overline{\tau }}^{\alpha \beta }\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}$$E9

The last term $$\overline{{\tau }^{\alpha \beta }\frac{\partial v{\text{'}}^{\beta }}{\partial x_{\alpha }}}=\overline{\rho }\epsilon$$of (6) is called « volumetric rate of dissipation for turbulence », due to the viscosity (and $$\epsilon$$is the rate of dissipation by unit of mass). This term is a transfer of energy from $$\tilde{k}$$to$$\tilde{e}$$, while the last term of (7) is the dissipation rate of kinetic energy of the mean motion, directly from $$\frac{1}{2}\overline{\rho }{\tilde{v}}^{\alpha }{\tilde{v}}^{\alpha }$$to$$\tilde{e}$$. There is in (6) a term corresponding to a transfer from$$\frac{1}{2}\overline{\rho }{\tilde{v}}^{\alpha }{\tilde{v}}^{\alpha }$$: the term $$-\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }}\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}$$is called the production term for$$\tilde{k}$$. The terms $$-\overline{p}\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\alpha }}-\overline{p\frac{\partial v{\text{'}}^{\alpha }}{\partial x_{\alpha }}}$$in (7) are also transfer terms from $$\frac{1}{2}\overline{\rho }{\tilde{v}}^{\alpha }{\tilde{v}}^{\alpha }$$and$$\tilde{k}$$, respectively, but these terms are reversible, unlike the last terms due to the viscosity. If the medium is incompressible (i.e.$$\frac{\partial v^{\alpha }}{\partial x_{\alpha }}\cong 0$$), these pressure terms are vanishing. When the velocity of the flow is low, very often the two last terms can be neglected in (7) with respect to the heat flux term (Chassaing et al., 2002).

Equations (1) to (3) are able to calculate the mean variables $$\overline{\rho },\tilde{e},{\tilde{v}}^{\alpha }$$if and only if the turbulent fluxes $$\overline{\rho },\tilde{e},{\tilde{v}}^{\alpha }$$and$$-\overline{p{v}^{\alpha }}+\overline{{\tau }^{\alpha \beta }{v}^{\beta }}-\overline{\rho v{\text{'}}^{\alpha }{e}_t}$$, additional to the classical viscous stress and heat conduction flux, can be provided by additional closure expressions, based on theoretical or empirical grounds. Instead of (3), we can use (7), and the transfer terms from $$\tilde{k}$$to $$\tilde{e}$$are also to be modelled. As we have previously said, this needs to take into account additional variables that represent in some sense the turbulent fluctuations, and, at least, $$\tilde{k}$$is one of these new variables. It is seen here that one balance equation for $$\tilde{k}$$can be found, which can be used provided again that the correlations involving the fluctuations that appear in this equation can be given also by expressions involving only$$\overline{\rho },\tilde{e},{\tilde{v}}^{\alpha }$$, $$\tilde{k}$$and eventually other ones…

The simplest model, the “Prandtl mixing length model”, does not use a balance equation for $$\tilde{k}$$but expresses it directly form the mean velocity gradient. We will see how that is done, and in addition it is postulated that turbulence involves a single (mean) turbulence length scale, $$\tilde{k}$$, which is simply related to the geometry of the turbulent flow. The second simplest turbulence model, the “Kolmogorov model”, uses the balance equation (6), the unknown terms in it being conveniently modelled with closure assumptions Then, the model provides expressions for the correlations involved in (1) to (6), in particular$$R_{\alpha \beta }$$, as functions of$$\overline{\rho },\tilde{e},{\tilde{v}}^{\alpha }$$, $$\tilde{k}$$, $$l_t$$and their spatial derivatives. The very popular “k-epsilon” model, instead of giving directly$$l_t$$, assume that the dissipation rate $$\epsilon$$itself can be calculated by an additional balance equation, postulates such an equation, and writes the correlations as functions of$$\overline{\rho },\tilde{e},{\tilde{v}}^{\alpha }$$, $$\tilde{k}$$, $$\epsilon$$and their spatial derivatives.

We will now show that the approach of irreversible thermodynamics can be of great help, once the type of model is chosen, to solve the problem of modelling. We will see that it justifies the usual choices done by researchers since 30 years, and suggests possibilities for unsolved problems.

2.2. Entropy for the Prandtl mixing length model

In this framework we expect to model the turbulent flow only with the variables$$\overline{\rho },\tilde{e},{\tilde{v}}^{\alpha }$$, given by their balance equations. It is assumed that the balance equation (6) for the kinetic energy of the fluctuating motion can be simplified in keeping only the production and destruction terms:

$$-\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }}\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}+\overline{p\frac{\partial v{\text{'}}^{\alpha }}{\partial x_{\alpha }}}=\overline{{\tau }^{\alpha \beta }\frac{\partial v{\text{'}}^{\beta }}{\partial x_{\alpha }}}$$E10

In addition, the rate of dissipation, which can be calculated through $$\tilde{k}$$and $$l_t$$only, is necessarily written as$$\overline{\rho }\epsilon =C_d{\tilde{k}}^{3/2}/l_t$$, where $$C_d$$is a non dimensional constant. This model has been initially proposed for incompressible flow, and the pressure term, with the divergence of velocity fluctuation, is totally negligible in this case, as well as the two pressure terms in (7). Then the turbulence kinetic energy is an algebraic function of the gradient of the mean velocity, and does not need to be calculated with the differential equation (6):

$$\overline{\rho }{C}_d\frac{{\tilde{k}}^{3/2}}{l_l}=-\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }}\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}$$E11

The three basic variables $$\overline{\rho },\tilde{e},{\tilde{v}}^{\alpha }$$being given, the entropy function has to be searched as $$\varsigma =\varsigma (\tilde{e},1/\overline{\rho })$$only (the velocity being a non-objective variable has not to appear), with the Gibbs and Euler-Gibbs equations classically written as:

$$\tilde{e}=\theta \varsigma -\pi \frac{1}{\overline{\rho }}+\mu$$,

$$d\tilde{e}=\theta d\varsigma -\pi d(\frac{1}{\overline{\rho }})$$E12

The two functions$$\theta$$, $$\pi$$are the partial derivatives of $$\tilde{e}$$with respect to $$\varsigma$$and$$1/\overline{\rho }$$, or $$1/\theta$$and $$\pi /\theta$$are the partial derivatives of $$\varsigma$$with respect to $$\tilde{e}$$and$$1/\overline{\rho }$$, and the relations $$1/\overline{\rho }$$et $$\pi =\pi (\tilde{e},1/\overline{\rho })$$are the laws of state, to be found. Of course, the formula $$\mu =\mu (\tilde{e},1/\overline{\rho })$$also can be deduced. In addition, these equations of state do satisfy necessarily the Maxwell relation:$$\mu =\mu (\tilde{e},1/\overline{\rho })$$.

The fluid itself, independently of the turbulent flow, has proper equations of state, and if we assume that it is an “ideal gas”, $$e=C_{v}T+e^0,p=\rho RT$$do hold, with three constants:$$C_v,e^0,R$$. Then, we can derive$$\tilde{e}=C_v\tilde{T}+e^0,\overline{p}=\overline{\rho }R\tilde{T}$$, just taking the average (weighted average for the first one, classic average for the second one). It appears that if we identify $$\theta =\tilde{T}$$and$$\pi =\overline{p}$$, we find$$\theta =\frac{\tilde{e}-e^0}{C_v},\pi =\overline{\rho }R\frac{\tilde{e}-e^0}{C_v}$$, and these relations are convenient equations of state because they satisfy the Maxwell relation : here$$\frac{\partial (1/\theta )}{\partial (1/\overline{\rho })}=\frac{\partial (\pi /\theta )}{\partial \tilde{e}}=0$$. The thermodynamic description of the model is complete. From these equations of state we can derive the entropy function, similarly as in the classical case of an ideal gas:$$\frac{\partial (1/\theta )}{\partial (1/\overline{\rho })}=\frac{\partial (\pi /\theta )}{\partial \tilde{e}}=0$$. We emphasize that this entropy is not the statistical mean value of the entropy of the fluid$$s=s(e,\rho )$$. Indeed, this mean value cannot be expressed in terms of $$\tilde{e},\overline{\rho }$$only, because $$s=s(e,\rho )$$is non linear, and is not usable in our modelling framework.

It has to be remarked that the previous development cannot be extended to the case of fluid satisfying a non linear Joule law$$e=e(T)$$. In this case, when the fluctuations are not infinitely small, it comes that$$\overline{\rho }\tilde{e}=f(\overline{\rho }\tilde{T},\overline{\rho T{\text{'}}^2},\overline{\rho T{\text{'}}^3},\mathrm{...})$$, and it is no more possible to limit the model to$$\varsigma =\varsigma (\tilde{e},1/\overline{\rho })$$. If the nonlinearity of $$e(T)$$is important, we have to consider, at least, the additional variable$$\overline{\rho T{\text{'}}^2}/\overline{\rho }$$.

Once we have defined the entropy and found the equations of state, we can use the classical approach of irreversible thermodynamics for deriving some laws for the additional fluxes appearing in the mean balance equations (2) and (7).

The Gibbs-Euler equation, now written as$$d\tilde{e}=\tilde{T}d\varsigma -\overline{p}d(\frac{1}{\overline{\rho }})$$, allows us to obtain the balance equation for the (new) entropy, showing the term of entropy production, which has to be always positive, and zero at equilibrium. Using (1), (2), (7), (8), one gets (with the “Lagrangian” derivative of the mean motion such that$$\frac{D}{Dt}(.)=\frac{\partial }{\partial t}(.)+{\tilde{v}}_{}^{\alpha }\frac{\partial }{\partial x_{\alpha }^{}}(.)$$):

$$\overline{\rho }\frac{D\varsigma }{Dt}=\frac{\partial }{\partial x_{\alpha }}(-\frac{{\overline{j}}_{\alpha }^Q+\overline{\rho v{\text{'}}^{\alpha }e}}{\tilde{T}})+({\overline{j}}_{\alpha }^Q+\overline{\rho v{\text{'}}^{\alpha }e})\frac{\partial 1/\tilde{T}}{\partial x_{\alpha }}+(\frac{{\overline{\tau }}_{}^{\alpha \beta }-\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }}}{\tilde{T}})\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}$$E13

The production terms are the two last terms, and in the last one the contribution due to the diagonal of the tensors is vanishing because the fluid is incompressible. It results that a simple extension of the linear classical irreversible thermodynamics is possible here, giving laws for the turbulent additional fluxes of heat and momentum.

Concerning the total friction term, a positive “total viscosity coefficient” $${\mu }_{tt}$$can be defined and:

$${\overline{\tau }}_{}^{\alpha \beta }-\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }}+\frac{2}{3}\overline{\rho }\tilde{k}{\delta }^{\alpha \beta }={\mu }_{tt}(\frac{\partial {\tilde{v}}_{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}_{\beta }}{\partial x_{\alpha }})$$E14

Considering that the fluid is Newtonian, one has already: $${\overline{\tau }}_{}^{\alpha \beta }=\mu (\frac{\partial {\tilde{v}}_{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}_{\beta }}{\partial x_{\alpha }})$$(here $$\mu =cst$$and$${\tilde{v}}^{\alpha }\cong {\overline{v}}^{\alpha }$$), and it appears for the pure turbulent contribution an “eddy viscosity coefficient” $${\mu }_t={\mu }_{tt}-\mu$$such as:

$$-\overline{\rho v{\text{'}}^{\alpha }v{\text{'}}^{\beta }}+\frac{2}{3}\overline{\rho }\tilde{k}{\delta }^{\alpha \beta }={\mu }_t(\frac{\partial {\tilde{v}}_{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}_{\beta }}{\partial x_{\alpha }})$$E15

This is the well known “Boussinesq relation” proposed for turbulent flows at the end of the 19th century.

Similarly, for the total heat flux, the linear law of irreversible thermodynamics gives:

$$-{\overline{j}}_{\alpha }^Q-\overline{\rho v{\text{'}}_{\alpha }e\text{'}}=(\lambda +{\lambda }_t)\frac{\partial \tilde{T}}{\partial x_{\alpha }}=(\lambda +\frac{{\mu }_t}{{\mathrm{Pr}}_t})\frac{\partial \tilde{T}}{\partial x_{\alpha }}$$E16

This introduces the classical « turbulent Prandtl number », constant, whose value (about 0.83) has been given by many experiments.

The eddy viscosity coefficient has now to be found. Just by dimensions, it can be written as$${\mu }_t=\overline{\rho }C{\tilde{k}}^{1/2}{l}_t$$, with (9) it comes:$$\tilde{k}=\frac{C}{C_d}{l}_t^2(\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }})\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}=\frac{C}{2C_d}{l}_t^2(\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }})(\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }})$$, and we get finally the well known formula of “Prandtl mixing length” :

$${\mu }_t=\overline{\rho }{C}_{ML}{l}_t^2{((\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }})(\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}))}^{1/2}$$E17

In this model, the turbulence length scale$$l_t$$, called “mixing length”, is given algebraically from the geometry of the flow. The constants $$C_{ML}$$and $$C_d$$have been found actual constants when the turbulent fluctuations are well developed, and then$${\mu }_t\gg \mu$$, but when it is not the case they vary in function of the “turbulence Reynolds number”$${\mathrm{Re}}_t={\mu }_t/\mu$$. Experiments have given such corrective functions.

The conclusion is then here that the classical approach of linear irreversible thermodynamics does justify perfectly the old empirical approach of Boussinesq and Prandtl.

2.3. Entropy for the K-epsilon model

We consider now the very popular “k-epsilon” model, where a balance equation for$$\tilde{k}$$, based on (6), and a balance equation for$${\epsilon }^{\ast }=\epsilon -2\nu \frac{\partial {\tilde{k}}^{1/2}}{\partial x_{\alpha }}\frac{\partial {\tilde{k}}^{1/2}}{\partial x_{\alpha }}$$, less firmly based on physics, are used. When turbulence is well established, $${\epsilon }^{\ast }\cong \epsilon$$, but the additional term is not negligible close to walls, where viscosity remains playing. There are other models that use equations for$$l_t$$, or$$1/{\tau }_t$$, the inverse of a time scale of turbulence, instead of$${\epsilon }^{\ast }$$, and similar discussions could be done for these models.

In this case, the independent variables $$\overline{\rho },\tilde{e},\tilde{k},{\epsilon }^{\ast }$$have to be taken into account in order to build an entropy function, and we have to search $$\varsigma =\varsigma (\tilde{e},1/\overline{\rho },\tilde{k},{\epsilon }^{\ast })$$with:

$$\tilde{e}=\theta \varsigma -\pi \frac{1}{\overline{\rho }}-y_k\tilde{k}-y_{\epsilon }{\epsilon }^{\ast }+\mu$$,

$$d\tilde{e}=\theta d\varsigma -\pi d(\frac{1}{\overline{\rho }})-y_{k}d\tilde{k}-y_{\epsilon }d{\epsilon }^{\ast }$$E18

Here$$\theta =\theta (\tilde{e},1/\overline{\rho },\tilde{k},{\epsilon }^{\ast })$$, $$\pi =\pi (\tilde{e},1/\overline{\rho },\tilde{k},{\epsilon }^{\ast })$$, $$y_k=y_k(\tilde{e},1/\overline{\rho },\tilde{k},{\epsilon }^{\ast })$$, $$y_{\epsilon }=y_{\epsilon }(\tilde{e},1/\overline{\rho },\tilde{k},{\epsilon }^{\ast })$$are the new equations of state to be found, satisfying the six corresponding Maxwell relations. We will have to discuss the finding of these equations of state in connexion with the discussion of the second principle. Concerning the $${\epsilon }^{\ast }$$equation, we use here the general following form, with a diffusion term and production and destruction terms:

$$\frac{\partial }{\partial t}(\overline{\rho }{\epsilon }^{\ast })+\frac{\partial }{\partial x_{\alpha }^{}}(\overline{\rho }{\tilde{v}}_{}^{\alpha }{\epsilon }^{\ast })=\frac{\partial }{\partial x_{\alpha }}(-J_{\epsilon }^{\alpha })-\overline{\rho }{C}_3\frac{2{\epsilon }^{\ast }}{3}\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\alpha }}-\overline{\rho }{C}_1\frac{{\epsilon }^{\ast }}{\tilde{k}}{R}_d^{\alpha \beta }\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}+C_4\frac{{\epsilon }^{\ast }}{\tilde{k}}\overline{p\frac{\partial v{\text{'}}^{\alpha }}{\partial x_{\alpha }}}-\overline{\rho }{C}_2\frac{\epsilon {\epsilon }^{\ast }}{\tilde{k}}$$E19

The C’s are constants (numbers, without dimension) or algebraic functions of the variables in some cases (the two first being always positive), depending of the exact version of the model used. Within the most classical k-ε model,$$C_1=1.44,C_2=1.92,C_3=1.44,C_4=0$$.

Using now (1), (6), (7) and (14), we can derive as in § 2.2 the balance equation for the entropy as:

$$\begin{array}{l}\overline{\rho }\frac{D\varsigma }{Dt}=\frac{\partial }{\partial x_{\alpha }}(-\frac{{\overline{j}}_{\alpha }^Q+\overline{\rho v{\text{'}}^{\alpha }e}}{\theta }-\frac{y_k{J}_k^{\alpha }+y_{\epsilon }{J}_{\epsilon }^{\alpha }}{\theta })+({\overline{j}}_{\alpha }^Q+\overline{\rho v{\text{'}}^{\alpha }e})\frac{\partial 1/\theta }{\partial x_{\alpha }}+\\ +J_k^{\alpha }\frac{\partial y_k/\theta }{\partial x_{\alpha }}+J_{\epsilon }^{\alpha }\frac{\partial y_{\epsilon }/\theta }{\partial x_{\alpha }}+\frac{1}{\theta }(\pi -\overline{p}-\frac{2}{3}\overline{\rho }\tilde{k}(y_k+\frac{C_3{\epsilon }^{\ast }{y}_{\epsilon }}{\tilde{k}}))\frac{\partial {\tilde{v}}^{\gamma }}{\partial x_{\gamma }}-\\ -\frac{1-y_k-C_4{\epsilon }^{\ast }{y}_{\epsilon }/\tilde{k}}{\theta }\overline{p\frac{\partial v{\text{'}}^{\beta }}{\partial x_{\beta }}}-\frac{y_k+C_1{\epsilon }^{\ast }{y}_{\epsilon }/\tilde{k}}{\theta }\overline{\rho }{R}_d^{\alpha \beta }\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}+\\ +\frac{1}{\theta }{\overline{\tau }}_d^{\alpha \beta }\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}+\frac{1}{3\theta }{\overline{\tau }}_{}^{\gamma \gamma }\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\alpha }}+\frac{1-y_k-C_2{y}_{\epsilon }{\epsilon }^{\ast }/\tilde{k}}{\theta }\overline{\rho }\epsilon \end{array}$$E20

We have written as $${\tau }_d^{\alpha \beta }$$and $$R_d^{\alpha \beta }$$the deviators of the tensors $${\tau }^{\alpha \beta }$$and$$R^{\alpha \beta }$$, and $$-J_k^{\alpha }$$is the diffusion term of (6). We recognize in this equation a diffusion term, the first on the right hand side, and all the following terms do constitute entropy production terms.

Before discussing the implications of the second principle, it is of interest to discuss the physical meaning of each term and to precise the equations of state. The first remark that has to be done is that the pressure has not to appear in the entropy production term, because it is not a source of irreversibility. Then the terms where the pressure appears have to be zero, and that implies that $$\pi =\overline{p}+\frac{2}{3}\overline{\rho }\tilde{k}{y}_k+\frac{2}{3}\overline{\rho }{C}_3{y}_{\epsilon }{\epsilon }^{\ast }$$and$$1-y_k-C_4\frac{{\epsilon }^{\ast }}{\tilde{k}}{y}_{\epsilon }=0$$.

In addition, within the last term of (15), the directly scalar part of the entropy production, there is one contribution of the energy equation but also the contribution $$-(y_k+C_2\frac{{\epsilon }^{\ast }}{\tilde{k}}{y}_{\epsilon })$$is due to the destruction terms of the equations for $$\tilde{k}$$and$${\epsilon }^{\ast }$$. Globally, this term has to be positive, but the contribution of $$\tilde{k}$$and $${\epsilon }^{\ast }$$may be dangerous for that. This problem is cancelled by choosing$$y_k+C_2\frac{{\epsilon }^{\ast }}{\tilde{k}}{y}_{\epsilon }=0$$, and that corresponds to a “mutual equilibrium of destructions rates” for $$\tilde{k}$$and$${\epsilon }^{\ast }$$. Assuming this mutual equilibrium ensures that no additional irreversibility is brought by this new kind of model.

These three relations do determine constraints on the equations of state, or on the C’s. It has to be remarked that the equation (14) has been built, and most often used, for cases where the turbulence is incompressible, then the divergence of the velocity fluctuations is vanishing, and the relation $$1-y_k-C_4\frac{{\epsilon }^{\ast }}{\tilde{k}}{y}_{\epsilon }=0$$does not matter in this circumstance.

Taking into account $$y_k+C_2\frac{{\epsilon }^{\ast }}{\tilde{k}}{y}_{\epsilon }=0$$and$$\pi =\overline{\rho }R\tilde{T}+\frac{2}{3}\overline{\rho }\tilde{k}{y}_k+\frac{2}{3}\overline{\rho }{C}_3{y}_{\epsilon }{\epsilon }^{\ast }$$, $$\tilde{e}=C_v\tilde{T}+e^0$$as well as the six Maxwell relations (saying that the second partial derivatives of $$\varsigma =\varsigma (\tilde{e},1/\overline{\rho },\tilde{k},{\epsilon }^{\ast })$$are identical when the derivation variables are commuted), it is possible with a few algebra to find a convenient set of four equations of state only in the case where $$C_3/C_2$$is a constant (possibly zero). The solution is:

$$\pi =\overline{\rho }R\theta$$, $$\tilde{e}-e^0=C_v\theta (1-\frac{2R_k}{3R}+\frac{2R_k{C}_3}{3R{C}_2})=C_v\tilde{T}$$, $$y_k=R_k\frac{\theta }{\tilde{k}}$$and

$$y_{\epsilon }=-\frac{R_k\theta }{C_2{\epsilon }^{\ast }}$$E21

,

where $$y_{\epsilon }=-\frac{R_k\theta }{C_2{\epsilon }^{\ast }}$$is a positive constant (with $$y_k$$chosen positive, without loss of generality). These formulas are not the unique solution of the problem, but they constitute a solution that can be supported by measurements, as approximations, at least in a certain limited range of variation of of the variables. One can remark that for keeping $$\theta$$positive, which is necessary in thermodynamics, it is necessary that$$1\frac{2R_k}{3R}(1-\frac{C_3}{C_2})$$.

We can now discuss in details the implications of the second principle, which prescribes that the entropy production rate, the right hand side of the following equation (16), has to be always positive, or zero at equilibrium.

$$\begin{array}{l}\overline{\rho }\frac{D\varsigma }{Dt}-\frac{\partial }{\partial x_{\alpha }}(\frac{{\overline{j}}_{\alpha }^Q+\overline{\rho v{\text{'}}^{\alpha }e}}{\theta }+\frac{y_k{J}_k^{\alpha }+y_{\epsilon }{J}_{\epsilon }^{\alpha }}{\theta })=({\overline{j}}_{\alpha }^Q+\overline{\rho v{\text{'}}^{\alpha }e})\frac{\partial 1/\theta }{\partial x_{\alpha }}+J_k^{\alpha }\frac{\partial y_k/\theta }{\partial x_{\alpha }}+J_{\epsilon }^{\alpha }\frac{\partial y_{\epsilon }/\theta }{\partial x_{\alpha }}\\ -\frac{1-y_k+y_k{C}_4/C_2}{\theta }\overline{p\frac{\partial v{\text{'}}^{\beta }}{\partial x_{\beta }}}-y_k\frac{1-C_1/C_2}{\theta }\overline{\rho }{R}_d^{\alpha \beta }\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}+\frac{1}{\theta }{\overline{\tau }}_d^{\alpha \beta }\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}+\frac{1}{3\theta }{\overline{\tau }}_{}^{\gamma \gamma }\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\alpha }}+\frac{1}{\theta }\overline{\rho }\epsilon \end{array}$$E22

-

We know already that $${\tau }_{}^{\gamma \gamma }={\mu }_v\frac{\partial v^{\alpha }}{\partial x_{\alpha }}$$and $${\tau }_d^{\alpha \beta }={\mu }_{}(\frac{\partial v^{\beta }}{\partial x_{\alpha }}+\frac{\partial v^{\beta }}{\partial x_{\alpha }}-\frac{{\delta }^{\alpha \beta }}{3}\frac{\partial v^{\gamma }}{\partial x_{\gamma }})$$with $${\mu }_v$$and $${\mu }_v$$positive, in such a way that $${\tau }^{\alpha \beta }\frac{\partial v^{\beta }}{\partial x_{\alpha }}$$is always positive, or zero if the velocity gradients are zero. We remember also that$$\overline{\rho }\epsilon =\overline{{\tau }^{\alpha \beta }\frac{\partial v{\text{'}}^{\beta }}{\partial x_{\alpha }}}$$, and then the three last term can be grouped giving$$\frac{1}{\theta }\overline{{\tau }^{\alpha \beta }\frac{\partial v^{\beta }}{\partial x_{\alpha }}}$$.Then, the group of the three last term of (16) is clearly always positive, by definition, even if the two first ones are not positive, depending on the fluctuations of density.

Secondly, the term $$-\frac{y_k}{\theta }(1-C_1/C_2)\overline{\rho }{R}_d^{\alpha \beta }\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}$$concerning the influence of turbulent friction has also to be positive always, or zero. At the same time, we remember that $$-\overline{\rho }{R}_d^{\alpha \beta }\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}$$is a production term for $$\tilde{k}$$and has to be essentially positive, although not necessarily everywhere. With $$C_2^{},C_1^{}$$constant, it is implied that $$1-C_1^{}/C_2^{}$$is necessarily positive and then the linear irreversible thermodynamics again supports the Boussinesq relation, here written for compressible flows:$$-\overline{\rho }{R}_d^{\alpha \beta }={\mu }_t(\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}+\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}-\frac{1}{3}\frac{\partial {\tilde{v}}^{\gamma }}{\partial x_{\gamma }}{\delta }^{\alpha \beta })$$, with $${\mu }_t$$positive.

Here again, we have to write$${\mu }_t=\overline{\rho }{C}_l{\tilde{k}}^{1/2}{l}_t=\overline{\rho }{C}_{\mu }{\tilde{k}}^2/\epsilon$$, $$C_{\mu }$$being positive, constant or function of the Reynolds number of turbulence, defined here as$${\mathrm{Re}}_t={\tilde{k}}^2/\epsilon \mu$$, if this number is sufficiently small. This is exactly what is used in all versions of the k-epsilon model. For versions able to consider small turbulence Reynolds number, $$1-C_1^{}/C_2^{}$$may become negative, and in this case this term has to be smaller than the three last terms corresponding to$$\frac{1}{\theta }\overline{{\tau }^{\alpha \beta }\frac{\partial v^{\beta }}{\partial x_{\alpha }}}\approx \frac{1}{\theta }\overline{{\tau }^{\alpha \beta }}\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }}$$.

The last point deals with the terms$$({\overline{j}}_{\alpha }^Q+\overline{\rho v{\text{'}}^{\alpha }e})\frac{\partial 1/\theta }{\partial x_{\alpha }}+J_k^{\alpha }\frac{\partial y_k/\theta }{\partial x_{\alpha }}+J_{\epsilon }^{\alpha }\frac{\partial y_{\epsilon }/\theta }{\partial x_{\alpha }}$$. Again, we see that the modelling of $$-J_k^{\alpha }$$and$$-J_{\epsilon }^{\alpha }$$, as well as $$-({\overline{j}}_{\alpha }^Q+\overline{\rho v{\text{'}}^{\alpha }e})$$can be done in terms of$$\frac{\partial \tilde{k}}{\partial x_{\alpha }}$$, $$\frac{\partial {\epsilon }^{\ast }}{\partial x_{\alpha }}$$and$$\frac{\partial \theta }{\partial x_{\alpha }}$$, with possible coupling. No coupling has been attempted so far, however. The coefficients of diffusivity of turbulence kinetic energy, dissipation rate, and heat are all proportional to $${\mu }_t$$previously defined.

As previously for the Prandtl mixing length model, we see that the well known k-epsilon model can be perfectly explained in the frame work of irreversible linear thermodynamics. The constraints we have found that $$C_2{C}_1$$do correspond with the common practice for large turbulence. The value of $$C_3$$is not well known, because the flows where $$\frac{\partial {\tilde{v}}_{\alpha }}{\partial x_{\alpha }}\ne 0$$are more complex, but we have shown that we need $$\frac{\partial {\tilde{v}}_{\alpha }}{\partial x_{\alpha }}\ne 0$$even if the C’s are varying. Considering$$C_4$$, we have found that$$1-R_k\frac{\theta }{\tilde{k}}(1-\frac{C_4}{C_2})=0$$. That implies then a constraint for$$C_4$$, which is$$C_4=C_2-\frac{\tilde{k}{C}_2}{\theta R_k}$$. In case where the influence of $$\tilde{k}/\theta$$is poor, $$C_4=C_2$$simply. However, it remains to find an approximation for $$\overline{p\frac{\partial v{\text{'}}_{\alpha }}{\partial x_{\alpha }}}$$when it is not negligible. In practice, this quantity plays only in flows where$$\frac{\partial v{\text{'}}_{\alpha }}{\partial x_{\alpha }}\ne 0$$, which are less well studied turbulent flows, and the experiments up to now have only confirmed that $$C_4$$has not to be put to zero...

3. Irreversible thermodynamics for two phase flows models

The theoretical description of two-phase flows is now firmly based by describing at each instant this medium as piece-wise continuous, and then considering that the model has to deal with statistical averages, like for turbulent flows (Drew, 1983). The randomness here is due both to the poorly known initial and limiting conditions concerning the presence of the phases, and to the existence of local perturbations permanently created by unstable phenomena within the flow, similarly to classical turbulent flows. There is again the need of defining how many mean variables are necessary, and of finding their equations. We first give here the set of variables and equations for representing the flow, without details concerning the calculations but with the relevant physical interpretation, and then we show how the Extended Irreversible Thermodynamics helps the modelling.

3.1. Variables and equations for two-phase flows modelling

The first classical attempt for representing a two-phase flow uses statistically averaged variables for volumetric mass, momentum, internal energy for the two phases, and their probability of presence, also known as the “volume fraction” of each phase. It is defined a “phase indicator” $${\Phi }_i(\overrightarrow{x},t)$$which is unity within the phase i and zero outside, and the statistical average $${\overline{\Phi }}_i(\overrightarrow{x},t)$$is nothing but the probability of presence of this phase. Of course,$${\displaystyle \sum _{i=1,2}}{\overline{\Phi }}_i(\overrightarrow{x},t)=1$$. For each phase, we define$${\overline{\rho }}_i=\overline{{\Phi }_i{\rho }_i}={\overline{\Phi }}_i{\overline{\rho }}_i^{\ast }$$, $$\overline{{\Phi }_i{\rho }_i{v}_i^{\alpha }}={\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{\tilde{v}}_i^{\alpha }$$,$$\overline{{\Phi }_i{\rho }_i{v}_i^{\alpha }}={\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{\tilde{v}}_i^{\alpha }$$.

The averaged forms of mass, momentum and energy equations can be found from the primitive equations for the piecewise continuous medium (Borghi, 2008). For momentum and energy they display additional diffusion fluxes due to the presence of local perturbations around the mean values, similarly as for the turbulent flows. We will call these terms “pseudo-turbulent” fluxes. The averaged equation for the mass of phase i is:

$$\frac{\partial }{\partial t}({\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i)+\frac{\partial }{\partial x_{\alpha }}({\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{\tilde{v}}_i^{\alpha })=-\overline{{\rho }_i{v}_{ri}^{\alpha }{n}_i^{\alpha }{\sigma }_s}$$E23

The last term on the R.H.S is the possible exchange of mass between the two phases: $$v_{ri}^{\alpha }$$is the velocity of the phase i relative to the interface, non zero only when there is exchange of mass between phases (vaporization-condensation for instance) ; $${\sigma }_s$$is the instantaneous interface area by unit of volume of the medium, mathematically defined by$$-n_i^{\alpha }{\sigma }_s=\frac{\partial {\Phi }_i}{\partial x_{\alpha }}$$, where $$n_i^{\alpha }$$is the normal to the interface, oriented outward phase i. Of course, the gradient of the phase indicator has the structure of a Dirac peak, and $${\sigma }_s$$is zero except on the interface (Kataoka, 1986). The averaged equation for the momentum of phase i is.

$$\frac{\partial }{\partial t}\left({\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{\tilde{v}}_i^{\beta }\right)+\frac{\partial }{\partial x_{\alpha }^{}}({\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{\tilde{v}}_i^{\alpha }{\tilde{v}}_i^{\beta })=\frac{\partial }{\partial x_{\alpha }}({\overline{\Phi }}_i{\overline{\sigma }}_i^{\ast \alpha \beta }-\overline{{\rho }_i{\Phi }_{i}v{\text{'}}_i^{\alpha }v{\text{'}}_i^{\beta }})+{\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{g}_{\beta }-\overline{{\rho }_i{v}_{ri}^{\alpha }{n}_i^{\alpha }{v}_i^{\beta }{\sigma }_s}+\overline{{\sigma }_i^{\alpha \beta }{n}_i^{\alpha }{\sigma }_s}$$E24

The pseudo turbulent momentum diffusion flux is$$-\overline{{\rho }_i{\Phi }_{}v{\text{'}}_i^{\alpha }v{\text{'}}_i^{\beta }}$$. The two last terms of (18) represent the momentum given to the phase i by the other one, by exchange of mass and by the contact force, respectively. The tensor $${\sigma }_i^{\alpha \beta }$$is the Cauchy stress tensor within the phase i. The equation for the mean total energy of phase i is:

$$\begin{array}{l}\frac{\partial }{\partial t}\left({\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{\tilde{e}}_{ti}\right)+\frac{\partial }{\partial x_{\alpha }^{}}({\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{\tilde{v}}_i^{\alpha }{\tilde{e}}_{ti})=\frac{\partial }{\partial x_{\alpha }}(-{\overline{\Phi }}_i{\overline{j_i^{Q\alpha }}}^{\ast }+{\overline{\Phi }}_i{\overline{{\sigma }_i^{\alpha \beta }}}^{\ast }{\tilde{v}}_i^{\beta }+\overline{{\Phi }_i{\sigma }_i^{\alpha \beta }v{\text{'}}_i^{\alpha }}-\overline{{\rho }_i{\Phi }_{i}v{\text{'}}_i^{\alpha }{e}_{ti}^{}})\\ +{\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{\tilde{v}}_i^{\beta }{g}_{\beta }-\overline{{\rho }_i{v}_{ri}^{\alpha }{n}_i^{\alpha }{e}_{ti}{\sigma }_s}+\overline{{\sigma }_i^{\alpha \beta }{n}_i^{\alpha }{v}_i^{\beta }{\sigma }_s}-\overline{j_i^{Q\alpha }{n}_i^{\alpha }{\sigma }_s}\end{array}$$E25

The total energy is defined as$$$$, the pseudo-turbulent energy flux is$${\tilde{e}}_{ti}={\tilde{e}}_i+\overline{{\rho }_i{\Phi }_i{v}_i^{\alpha }{v}_i^{\alpha }}/2{\overline{\rho }}_i={\tilde{e}}_i+{\tilde{v}}_i^{\alpha }{\tilde{v}}_i^{\alpha }/2+\overline{{\rho }_i{\Phi }_{i}v{\text{'}}_i^{\alpha }v{\text{'}}_i^{\alpha }}/2{\overline{\rho }}_i$$, and the last three terms are the exchanges of energy from the other phase, due to exchange of mass, power of contact force, and heat flux, respectively.

It is possible also to write a balance equation for the volume of phases, simply obtained from the convection equation of the field of$$-\overline{{\rho }_i{\Phi }_{}v{\text{'}}_i^{\alpha }{e}_{ti}}$$:$${\Phi }_i(\overrightarrow{x},t)$$. When averaged, this equation gives:

$$\frac{\partial }{\partial t}{\Phi }_i+v_s^{\alpha }\frac{\partial }{\partial x_{\alpha }}{\Phi }_i=\frac{\partial }{\partial t}{\Phi }_i+(v_i^{\alpha }-v_{ri}^{\alpha })\frac{\partial }{\partial x_{\alpha }}{\Phi }_i=0$$E27

The terms on the R.H.S. have to be approximated by models, the second one being only due to the exchange of mass between phases.

There are also important instantaneous interface conservation relations, linking at the interface the exchanges between the phases: first, the mass lost or gained by one phase is identical to the mass gained or lost by the other phase; second, if we neglect the surface tension phenomena, it follows that the momentum lost or gained by one phase is identical to the momentum gained or lost by the other one ; and third a similar relation holds for the total energy (Kataoka, 1986). That leads to:

$$\frac{\partial }{\partial t}\left({\overline{\Phi }}_i\right)+{\tilde{v}}_i^{\alpha }\frac{\partial }{\partial x_{\alpha }^{}}({\overline{\Phi }}_i)=-\overline{v_i^{\text{'}\alpha }\frac{\partial {\Phi }_i}{\partial x_{\alpha }}}+\overline{v_{ri}^{\alpha }\frac{\partial {\Phi }_i}{\partial x_{\alpha }}}$$E28

At the interface, the tangential velocity and the temperature fields are continuous, but there are discontinuities of density, momentum and energy per unit mass. The equations (21) give “jump relations” relating the values of these quantities on both sides of the discontinuity. We will not consider further here the exchanges of mass between phases, and these relations are simplified as:

$${\displaystyle \sum _{i=1,2}}{({\rho }_i{v}_{ri}^{\alpha }{n}_i^{\alpha })}_s=0,{\displaystyle \sum _{i=1,2}}{(-{\rho }_i{v}_{ri}^{\alpha }{n}_i^{\alpha }{v}_i^{\beta }+{\sigma }_i^{\alpha \beta }{n}_i^{\alpha })}_s=0,{\displaystyle \sum _{i=1,2}}{(-{\rho }_i{v}_{ri}^{\alpha }{n}_i^{\alpha }{e}_{ti}+{\sigma }_i^{\alpha \beta }{n}_i^{\alpha }{v}_i^{\beta }-j_i^{Q\alpha }{n}_i^{\alpha })}_s=0$$E29

The third relation of (22) is deduced from the second one because both phases have the same velocity on the interface.

The instantaneous behaviour of each phase has to be given by the knowledge of $$v_{r1}^{\alpha }=v_{r2}^{\alpha }=0,{({\sigma }_1^{\alpha \beta }{n}_1^{\alpha })}_s+{({\sigma }_2^{\alpha \beta }{n}_2^{\alpha })}_s=0,{({\sigma }_1^{\alpha \beta }{n}_1^{\alpha }{v}_1^{\beta })}_s+{({\sigma }_2^{\alpha \beta }{n}_2^{\alpha }{v}_2^{\beta })}_s=0,{(j_1^{Q\alpha }{n}_1^{\alpha })}_s+{(j_2^{Q\alpha }{n}_2^{\alpha })}_s=0$$and $${\sigma }_i^{\alpha \beta }$$as function of$$j_i^{Q\alpha }$$. We will notice here $${\rho }_i,v_i^{\alpha },T_i$$for any type of phase (but usually for a solid phase the pressure is defined by the third of the trace of$${\sigma }_i^{\alpha \beta }=-p_i{\delta }^{\alpha \beta }+{\tau }_i^{\alpha \beta }$$, with an opposite sign). For a fluid phase, $${\sigma }_i^{\alpha \beta }$$is the viscous stress tensor, which can be assumed Newtonian for instance, and an equation of state gives$${\tau }_i^{\alpha \beta }$$. For solid particles, we are not concerned with the full elasticity properties and we can limit the description as a “quasi-liquid”, very viscous. For any type of phase the Fourier law for heat conduction can be adopted and a Joule law linking the temperature to the internal energy can be considered, neglecting the energetic aspects of compressibility for liquid or solid phases. Very often it suffices to take$$p_i=p_i({\rho }_i,T_i)$$.

One sees the similarities of the problem of modelling these equations with the one of modelling turbulence in single phase flows. We will show now how the irreversible thermodynamics can help this modelling task. For simplicity, we will consider in the sequel that there is no exchange of mass between the phases, this could be taken into account later.

3.2. Entropy for a simple two-phase flows model

We will consider the simple case of a model similar to the Prandtl mixing length model, where simply the variables: $${\rho }_i=cst$$have to be represented by their own balance equations. The development then follows the same path as § 2.2 above, but with two phases i=1, i=2.

The kinetic energies of fluctuations, namely$${\overline{\rho }}_i^{\ast },{\tilde{v}}_i^{\alpha },{\tilde{e}}_i,{\overline{\Phi }}_i,i=1,2$$, follow a balance equation that can be found from the primitive equations. It writes, without exchange of mass:

$$k_i=+\overline{{\rho }_i{\Phi }_{i}v{\text{'}}_i^{\alpha }v{\text{'}}_i^{\alpha }}/2{\overline{\rho }}_i,i=1,2$$E30

We do not consider here, in the Prandtl mixing length framework, that each $$\begin{array}{l}\frac{\partial }{\partial t}({\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{k}_i)+\frac{\partial }{\partial x_{\alpha }}({\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{k}_i{\tilde{v}}_i^{\alpha })=\frac{\partial }{\partial x_{\alpha }}(-\overline{{\rho }_i{\Phi }_{i}v{\text{'}}_i^{\beta }v{\text{'}}_i^{\beta }v{\text{'}}_{}^{\alpha }/2}+\overline{{\Phi }_{i}v{\text{'}}_i^{\beta }{\sigma }_i^{\alpha \beta }})-\overline{{\rho }_i{\Phi }_{i}v{\text{'}}_i^{\alpha }v{\text{'}}_i^{\beta }}\frac{\partial {\tilde{v}}_i^{\beta }}{\partial x_{\alpha }}\\ +\overline{p_i{\Phi }_i\frac{\partial v{\text{'}}_i^{\alpha }}{\partial x_{\alpha }}}-{\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{\epsilon }_i+\overline{{\sigma }_i^{\alpha \beta }{n}_i^{\beta }v{\text{'}}_i^{\alpha }{\sigma }_s}\end{array}$$does follow this full equation, we assume that it can be reduced to a local balance between production terms and destructions terms, giving simply that :

$$k_i$$E31

The last term of (23) and (23’) is due to the exchange of energy from the other phase, by the contact force.

With the point of view of the Prandtl model applied to both phases, we have to build an entropy with the variables$$-\overline{{\rho }_i{\Phi }_{i}v{\text{'}}_i^{\alpha }v{\text{'}}_i^{\beta }}\frac{\partial {\tilde{v}}_i^{\beta }}{\partial x_{\alpha }}+\overline{p_i{\Phi }_i\frac{\partial v{\text{'}}_i^{\alpha }}{\partial x_{\alpha }}}+\overline{{\sigma }_i^{\alpha \beta }{n}_i^{\beta }v{\text{'}}_i^{\alpha }{\sigma }_s}={\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{\epsilon }_i$$, such as$${\overline{\rho }}_i^{\ast },{\tilde{v}}_i^{\alpha },{\tilde{e}}_i,{\overline{\Phi }}_i,i=1,2$$, because entropy is an extensive quantity. The entropies of each phase are functions only of the variables of this phase, unlike molecular mixtures:$$\overline{\rho }\varsigma ={\displaystyle \sum _{i=1,2}{\overline{\rho }}_i^{\ast }}{\overline{\Phi }}_i{\varsigma }_i$$. As in § 2.2, we can adopt

$${\varsigma }_i={\varsigma }_i({\overline{\rho }}_i^{\ast },{\tilde{e}}_i)$$,

$${\tilde{e}}_i={\tilde{T}}_i{\varsigma }_i-{\overline{p}}_i^{\ast }\frac{1}{{\overline{\rho }}_i^{\ast }}+{\mu }_i$$E32

Similarly, we take linear equations of state for phases, which can be averaged without needing additional variables. For a gaseous phases$$d{\tilde{e}}_i={\tilde{T}}_{i}d{\varsigma }_i-{\overline{p}}_i^{\ast }d(\frac{1}{{\overline{\rho }}_i^{\ast }})$$, and for liquid or solid phases$${\tilde{e}}_i=C_{vi}{\tilde{T}}_i+e_i^0,{\overline{p}}_i^{\ast }={\overline{\rho }}_i^{\ast }{R}_i{\tilde{T}}_i$$, neglecting the influence of compressibility on the internal energy.

We can first show that the convective time derivative for the mean entropy of the medium is such that:

$${\tilde{e}}_s=C_s{\tilde{T}}_s+e_s^0,{\overline{p}}_s^{\ast }=K_s({\overline{\rho }}_s^{\ast }-{\rho }_0)$$E33

The convective time derivative for the medium follows the averaged velocity, while for the phases they follow the averaged velocity of each phase. Of course, the Euler-Gibbs relation is written now as:

$$\overline{\rho }\frac{D\varsigma }{Dt}={\displaystyle \sum _i{\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i\frac{D{\varsigma }_i}{D{t}_i}}-\frac{\partial }{\partial x_{\alpha }}({\displaystyle \sum _i{\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i}{\varsigma }_i({\tilde{v}}_i^{\alpha }-{\tilde{v}}^{\alpha }))$$E34

It is easy to verify that the equation (17) allows to write:

$$\frac{D}{D{t}_i}{\tilde{e}}_i={\tilde{T}}_i\frac{D}{D{t}_i}{\varsigma }_i-{\overline{p}}_i^{\ast }\frac{D}{D{t}_i}(\frac{1}{{\overline{\rho }}_i^{\ast }})$$, without exchange of mass between phases. By using (20), we obtain as well:

$${\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i\frac{D}{D{t}_i}(\frac{1}{{\overline{\rho }}_i^{\ast }})=-\frac{{\overline{\Phi }}_i}{{\overline{\rho }}_i^{\ast }}\frac{D{\overline{\rho }}_i^{\ast }}{D{t}_i}={\overline{\Phi }}_i\frac{\partial {\tilde{v}}_i^{\alpha }}{\partial x_{\alpha }}+\frac{D{\overline{\Phi }}_i}{D{t}_i}$$E35

It is necessary then to obtain the balance equation for the mean internal energy of phases from (19), removing the kinetic energy of the mean motion by taking into account (18), and removing also the kinetic energy of the fluctuations with (23). It is found:

$${\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i\frac{D}{D{t}_i}(\frac{1}{{\overline{\rho }}_i^{\ast }})={\overline{\Phi }}_i\frac{\partial {\tilde{v}}_i^{\alpha }}{\partial x_{\alpha }}-\overline{v{\text{'}}_i^{\alpha }\frac{\partial {\Phi }_i}{\partial x_{\alpha }}}$$E36

With (23’), it gives:

$$\begin{array}{l}{\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i\frac{D{\tilde{e}}_i}{D{t}_i}=\frac{\partial }{\partial x_{\alpha }}(-{\overline{\Phi }}_i{\overline{j_i^{Q\alpha }}}^{\ast }-\overline{{\rho }_i{\Phi }_i{e}_{i}v{\text{'}}_i^{\alpha }})-{\overline{\Phi }}_i{\overline{p}}_i^{\ast }\frac{\partial {\tilde{v}}_i^{\alpha }}{\partial x_{\alpha }}\\ +{\overline{\Phi }}_i{\overline{{\tau }_i^{\alpha \beta }}}^{\ast }\frac{\partial {\tilde{v}}_i^{\beta }}{\partial x_{\alpha }}-\overline{p_i{\Phi }_i\frac{\partial v{\text{'}}_i^{\beta }}{\partial x_{\beta }}}-\overline{j_i^{Q\alpha }{n}_i^{\alpha }{\sigma }_s}+{\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i{\epsilon }_i\end{array}$$

$${\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i\frac{D{\tilde{e}}_i}{D{t}_i}=\frac{\partial }{\partial x_{\alpha }}(-{\overline{\Phi }}_i{\overline{j_i^{Q\alpha }}}^{\ast }-\overline{{\rho }_i{\Phi }_i{e}_{i}v{\text{'}}_i^{\alpha }})-{\overline{\Phi }}_i{\overline{p}}_i^{\ast }\frac{\partial {\tilde{v}}_i^{\alpha }}{\partial x_{\alpha }}+{\overline{\Phi }}_i({\overline{{\tau }_i^{\alpha \beta }}}^{\ast }-{\overline{\rho }}_i^{\ast }{R}_i^{\alpha \beta })\frac{\partial {\tilde{v}}_i^{\beta }}{\partial x_{\alpha }}$$E37

,

where$${\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i\frac{D{\tilde{e}}_i}{D{t}_i}=\frac{\partial }{\partial x_{\alpha }}(-{\overline{\Phi }}_i{\overline{j_i^{Q\alpha }}}^{\ast }-\overline{{\rho }_i{\Phi }_i{e}_{i}v{\text{'}}_i^{\alpha }})-{\overline{\Phi }}_i{\overline{p}}_i^{\ast }\frac{\partial {\tilde{v}}_i^{\alpha }}{\partial x_{\alpha }}+{\overline{\Phi }}_i({\overline{{\tau }_i^{\alpha \beta }}}^{\ast }-{\overline{\rho }}_i^{\ast }{R}_i^{\alpha \beta })\frac{\partial {\tilde{v}}_i^{\beta }}{\partial x_{\alpha }}$$is the deviator of the Reynolds tensor for phase i. The two last terms take into account the exchange of heat by conduction $$R_{id}^{\alpha \beta }$$and the exchange of kinetic energy by contact force$${\overline{\phi }}_i^Q=-\overline{j_i^{Q\alpha }{n}_i^{\alpha }{\sigma }_s}$$, respectively, from the other phase. In the momentum balance (18) the total contact force has to be split in$$\overline{{\sigma }_i^{\alpha \beta }{n}_i^{\beta }v{\text{'}}_i^{\alpha }{\sigma }_s}$$. The first part is simply the influence of the variation of the volume of the phase, while $$\overline{{\sigma }_i^{\alpha \beta }{n}_i^{\beta }{\sigma }_s}=+\overline{p}\frac{\partial {\overline{\Phi }}_i}{\partial x_{\alpha }}+{\overline{F}}_{si}^{\times \beta }$$is really the force from the other phase, and the jump condition (22) gives$${\overline{F}}_{si}^{\times \beta }=\overline{{\tau }_i^{\times \alpha \beta }{n}_i^{\alpha }{\sigma }_s}$$. It follows that$${\displaystyle \sum _{i=1,2}}{\overline{F}}_{si}^{\times \beta }=0$$, giving finally:

$$\overline{{\sigma }_i^{\alpha \beta }{n}_i^{\beta }v{\text{'}}_i^{\alpha }{\sigma }_s}=+\overline{p}\overline{v{\text{'}}_i^{\alpha }\frac{\partial {\Phi }_i}{\partial x_{\alpha }}}+\overline{F_{si}^{\times \alpha }v{\text{'}}_i^{\alpha }}$$E38

Now from (24), (25), (26), (27), we obtain the entropy balance equation for the medium:

$$\overline{{\sigma }_i^{\alpha \beta }{n}_i^{\beta }v{\text{'}}_i^{\alpha }{\sigma }_s}=+\overline{p}\overline{v{\text{'}}_i^{\alpha }\frac{\partial {\Phi }_i}{\partial x_{\alpha }}}+\overline{F_{si}^{\times \alpha }v{\text{'}}_i^{\alpha }}$$E39

We recognize in this equation the rate of production of entropy, which has to be positive or zero in any case.

The last sum deals with the exchange of heat between the phases. As$$\begin{array}{l}\overline{\rho }\frac{D\varsigma }{Dt}=\frac{\partial }{\partial x_{\alpha }}(-{\displaystyle \sum _i({\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i}{\varsigma }_i({\tilde{v}}_i^{\alpha }-{\tilde{v}}^{\alpha })+\frac{{\overline{\Phi }}_i}{{\tilde{T}}_i}{\overline{j_i^{Q\alpha }}}^{\ast }+\frac{\overline{{\rho }_i{\Phi }_i{e}_{i}v{\text{'}}_i^{\alpha }}}{{\tilde{T}}_i}))+{\displaystyle \sum _i-}({\overline{\Phi }}_i{\overline{j_i^{Q\alpha }}}^{\ast }+\overline{{\rho }_i{\Phi }_i{e}_{i}v{\text{'}}_i^{\alpha }})\\ \frac{\partial }{\partial x_{\alpha }}(\frac{1}{{\tilde{T}}_i})+{\displaystyle \sum _i\frac{{\overline{\Phi }}_i}{{\tilde{T}}_i}}({\overline{{\tau }_i^{\alpha \beta }}}^{\ast }-{\overline{\rho }}_i^{\ast }{R}_i^{\alpha \beta })\frac{\partial {\tilde{v}}_i^{\beta }}{\partial x_{\alpha }}+{\displaystyle \sum _i\frac{\overline{F_{si}^{\times \beta }v{\text{'}}_i^{\beta }}}{{\tilde{T}}_i}}-{\displaystyle \sum _i\frac{{\overline{p}}_i^{\ast }-\overline{p}}{{\tilde{T}}_i}}\overline{v{\text{'}}_i^{\alpha }\frac{\partial {\Phi }_i}{\partial x_{\alpha }}}+{\displaystyle \sum _i\frac{{\overline{\phi }}_i^Q}{{\tilde{T}}_i}}\end{array}$$, this term can be written as:$$\begin{array}{l}\overline{\rho }\frac{D\varsigma }{Dt}=\frac{\partial }{\partial x_{\alpha }}(-{\displaystyle \sum _i({\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i}{\varsigma }_i({\tilde{v}}_i^{\alpha }-{\tilde{v}}^{\alpha })+\frac{{\overline{\Phi }}_i}{{\tilde{T}}_i}{\overline{j_i^{Q\alpha }}}^{\ast }+\frac{\overline{{\rho }_i{\Phi }_i{e}_{i}v{\text{'}}_i^{\alpha }}}{{\tilde{T}}_i}))+{\displaystyle \sum _i-}({\overline{\Phi }}_i{\overline{j_i^{Q\alpha }}}^{\ast }+\overline{{\rho }_i{\Phi }_i{e}_{i}v{\text{'}}_i^{\alpha }})\\ \frac{\partial }{\partial x_{\alpha }}(\frac{1}{{\tilde{T}}_i})+{\displaystyle \sum _i\frac{{\overline{\Phi }}_i}{{\tilde{T}}_i}}({\overline{{\tau }_i^{\alpha \beta }}}^{\ast }-{\overline{\rho }}_i^{\ast }{R}_i^{\alpha \beta })\frac{\partial {\tilde{v}}_i^{\beta }}{\partial x_{\alpha }}+{\displaystyle \sum _i\frac{\overline{F_{si}^{\times \beta }v{\text{'}}_i^{\beta }}}{{\tilde{T}}_i}}-{\displaystyle \sum _i\frac{{\overline{p}}_i^{\ast }-\overline{p}}{{\tilde{T}}_i}}\overline{v{\text{'}}_i^{\alpha }\frac{\partial {\Phi }_i}{\partial x_{\alpha }}}+{\displaystyle \sum _i\frac{{\overline{\phi }}_i^Q}{{\tilde{T}}_i}}\end{array}$$, and the linear classical irreversible thermodynamics proposes$$\begin{array}{l}\overline{\rho }\frac{D\varsigma }{Dt}=\frac{\partial }{\partial x_{\alpha }}(-{\displaystyle \sum _i({\overline{\rho }}_i^{\ast }{\overline{\Phi }}_i}{\varsigma }_i({\tilde{v}}_i^{\alpha }-{\tilde{v}}^{\alpha })+\frac{{\overline{\Phi }}_i}{{\tilde{T}}_i}{\overline{j_i^{Q\alpha }}}^{\ast }+\frac{\overline{{\rho }_i{\Phi }_i{e}_{i}v{\text{'}}_i^{\alpha }}}{{\tilde{T}}_i}))+{\displaystyle \sum _i-}({\overline{\Phi }}_i{\overline{j_i^{Q\alpha }}}^{\ast }+\overline{{\rho }_i{\Phi }_i{e}_{i}v{\text{'}}_i^{\alpha }})\\ \frac{\partial }{\partial x_{\alpha }}(\frac{1}{{\tilde{T}}_i})+{\displaystyle \sum _i\frac{{\overline{\Phi }}_i}{{\tilde{T}}_i}}({\overline{{\tau }_i^{\alpha \beta }}}^{\ast }-{\overline{\rho }}_i^{\ast }{R}_i^{\alpha \beta })\frac{\partial {\tilde{v}}_i^{\beta }}{\partial x_{\alpha }}+{\displaystyle \sum _i\frac{\overline{F_{si}^{\times \beta }v{\text{'}}_i^{\beta }}}{{\tilde{T}}_i}}-{\displaystyle \sum _i\frac{{\overline{p}}_i^{\ast }-\overline{p}}{{\tilde{T}}_i}}\overline{v{\text{'}}_i^{\alpha }\frac{\partial {\Phi }_i}{\partial x_{\alpha }}}+{\displaystyle \sum _i\frac{{\overline{\phi }}_i^Q}{{\tilde{T}}_i}}\end{array}$$.

The second last term involves$${\overline{\phi }}_1^Q\propto (\frac{1}{{\tilde{T}}_1}-\frac{1}{{\tilde{T}}_2})$$, expression implied in the closure of (20), and we see that assuming that $$\overline{v{\text{'}}_i^{\alpha }\frac{\partial {\Phi }_i}{\partial x_{\alpha }}}$$eliminates this term and at the same time the need of (20). Instead of this assumption, the linear irreversible thermodynamics suggests the hypotheses that $${\overline{p}}_i^{\ast }=\overline{p}$$and gives a closure for (20). The phase volume fraction can be calculated then with$$-\overline{v{\text{'}}_i^{\alpha }\frac{\partial {\Phi }_i}{\partial x_{\alpha }}}\propto ({\overline{p}}_i^{\ast }-\overline{p})$$. Such an equation is almost identical to the equation proposed by M. Baer and I. Nunziato for compressed granular medium (Baer & Nunziato, 1986).

The two first terms of the entropy production deal with the diffusion fluxes of heat and momentum. For momentum, when each phase is divergence free (which is originally the scope of Prandtl model), we find for each phase the same formulation as for a single phase turbulent flow in §2.2. Similarly, for heat diffusion, if this term is considered uncoupled with others, the classical model found in §2.2 is again obtained for each phase.

The term $$\frac{\partial }{\partial t}\left({\overline{\Phi }}_i\right)+{\tilde{v}}_i^{\alpha }\frac{\partial }{\partial x_{\alpha }^{}}({\overline{\Phi }}_i)=-\overline{v{\text{'}}_i^{\alpha }\frac{\partial {\Phi }_i}{\partial x_{\alpha }}}=K({\overline{p}}_i^{\ast }-\overline{p})$$is connected with the exchange of momentum between phases. The “jump relations” (22), similarly as$${\displaystyle \sum _i\frac{\overline{F_{si}^{\times \beta }v{\text{'}}_i^{\beta }}}{{\tilde{T}}_i}}={\displaystyle \sum _i\frac{1}{{\tilde{T}}_i}}(\overline{F_{si}^{\times \beta }{v}_i^{\beta }}-{\overline{F}}_{si}^{\times \beta }{\tilde{v}}_i^{\beta })$$, gives $${\displaystyle \sum _{i=1,2}}{\overline{F}}_{si}^{\times \beta }=0$$(and gave$${\displaystyle \sum _{i=1,2}}\overline{F_{si}^{\times \beta }{v}_i^{\beta }}=0$$previously used). Because at the interface between the phases the velocities of phases are equal, an averaged interface velocity $${\displaystyle \sum _{i=1,2}{\overline{\phi }}_i^Q}=0$$can be defined such that$${\overline{v}}_s^{\beta }$$, and the related source term of entropy can be written$$\overline{F_{si}^{\times \beta }{v}_i^{\beta }}={\overline{F}}_{si}^{\times \beta }{\overline{v}}_s^{\beta }$$. The classical linear irreversible thermodynamics leads to$${\displaystyle \sum _i\frac{1}{{\tilde{T}}_i}}{\overline{F}}_{si}^{\times \beta }({\overline{v}}_s^{\beta }-{\tilde{v}}_i^{\beta })$$, and with$${\displaystyle \sum _i\frac{1}{{\tilde{T}}_i}}{\overline{F}}_{si}^{\times \beta }({\overline{v}}_s^{\beta }-{\tilde{v}}_i^{\beta })$$, that gives an expression for the averaged interface velocity :$${\overline{F}}_{s1}^{\times \beta }=-{\overline{F}}_{s2}^{\times \beta }$$. The$${\overline{v}}_s^{\beta }=\frac{K_1{\tilde{v}}_1^{\beta }+K_2{\tilde{v}}_2^{\beta }}{K_1+K_2}$$are positive constants, depending on the characteristics of the two phases in contact, to be found with experiments. If we consider the possibility of a coupling, in agreement with the so called Curie theorem, we find that it is possible that the contact force $$K_i$$to the phase i is linked to the gradients of temperature within the phases. We recover here a “thermophoretic “effect.

All these discussions do justify the assumptions currently done for the modelling of two-phase flows, in a large range of applications.

4. Irreversible thermodynamics for modelling of granular media

Granular flows are nothing but two-phase flows of solid irregular particles within a gaseous or liquid phase. But the high number of particles per unit of volume, and the low velocity that may persist allow numerous and long duration contacts between the particles, unlikely to usual two-phase flows. The approach that has been presented in §3 for two phase flows can be used also, but the problem that arises immediately is that very often the solid phase, locally, is built with two or three particles (or more) in contact, and consequently the solid here is not a true solid for which behaviour laws are known. Actually, each grain is one solid, but a couple of grains does not. So, it is necessary to consider the medium as a multi-phase medium, where each grain represents one phase. This approach has been proposed recently by R.Borghi and S.Bonelli (Borghi § Bonelli, 2007), and we will develop a while its connection with Irreversible Thermodynamics.

4.1. Basic equations for granular media

For granular media, the equations (17) to (20) do hold, but we can simplify them because it is not needed to follow the motion of each grain. Considering that there are N solid grains, there are N+1 phases and N+1$${\overline{F}}_{si}^{\times \beta }$$, but we are interested only in the knowledge of$${\overline{\rho }}_i^{\ast },{\tilde{v}}_i^{\alpha }$$, and$${\overline{\Phi }}_g={\displaystyle \sum _{i=1,N}}{\overline{\Phi }}_i,{\overline{\rho }}_g^{\ast }{\overline{\Phi }}_g={\displaystyle \sum _{i=1,N}}\overline{{\Phi }_i{\rho }_i},{\overline{\rho }}_g^{\ast }{\tilde{v}}_g^{\alpha }{\overline{\Phi }}_g={\displaystyle \sum _{i=1,N}}\overline{{\rho }_i{\Phi }_i{v}_i^{\alpha }}$$. From (17), (18), we get the mass and momentum balance equations of the “grains phase” as:

$${\overline{\rho }}_f^{\ast }\equiv {\overline{\rho }}_{N+1}^{\ast },{\tilde{v}}_f^{\alpha }\equiv {\tilde{v}}_{N+1}^{\alpha }$$E40

Here the mean intrinsic density of the grains $$\frac{\partial }{\partial t}({\overline{\rho }}_g^{\ast }{\overline{\Phi }}_g)+\frac{\partial }{\partial x_{\alpha }}({\overline{\rho }}_g^{\ast }{\overline{\Phi }}_g{\tilde{v}}_g^{\alpha })=0$$can be considered generally as known and constant. If all the grains are made with the same almost incompressible solid, $${\overline{\rho }}_g^{\ast }$$is the density of this solid$${\overline{\rho }}_g^{\ast }$$.

$${\overline{\rho }}_g^{\ast }={\rho }_s$$E41

We see that the effective Cauchy stress tensor for the grains has two parts:$$\begin{array}{l}\frac{\partial }{\partial t}\left({\overline{\rho }}_g^{\ast }{\overline{\Phi }}_g{\tilde{v}}_g^{\beta }\right)+\frac{\partial }{\partial x_{\alpha }^{}}({\overline{\rho }}_g^{\ast }{\overline{\Phi }}_g{\tilde{v}}_g^{\alpha }{\tilde{v}}_g^{\beta })=\frac{\partial }{\partial x_{\alpha }}({\displaystyle \sum _{i=1,N}}{\overline{\Phi }}_i{\overline{\sigma }}_i^{\ast \alpha \beta }-{\displaystyle \sum _{i=1,N}}\overline{{\rho }_i{\Phi }_{i}v\text{'}{\text{'}}_i^{\alpha }v\text{'}{\text{'}}_i^{\beta }})+{\overline{\rho }}_g^{\ast }{\overline{\Phi }}_g{g}_{\beta }\\ +{\displaystyle \sum _{i=1,N}}\overline{{\sigma }_i^{\alpha \beta }{n}_i^{\alpha }{\sigma }_s^i}\end{array}$$. The first one takes into account the stresses within the grains but mainly also the influences of the lasting contacts, with friction between grains depending on pertinent laws, and we can call it “mean Cauchy stress of lasting contacts”; the second one is due to the velocity fluctuations, here defined with respect to the averaged velocity of all the grains:$${\overline{\sigma }}_{tg}^{\alpha \beta }={\displaystyle \sum _{i=1,N}}{\overline{\Phi }}_i{\overline{\sigma }}_i^{\ast \alpha \beta }-{\displaystyle \sum _{i=1,N}}\overline{{\rho }_i{\Phi }_{i}v\text{'}{\text{'}}_i^{\alpha }v\text{'}{\text{'}}_i^{\beta }}$$, and is similar to the Reynolds tensor for turbulent flows, and we will call it “mean kinetic stress”. It is important only when the motion of the medium is well established, but even with slow mean motion there are some large and sudden fluctuations of velocity. The modelling of both terms has to be studied with care. Without models, eq.(31) is without interest.

The last term of (31) is the influence of the fluid phase. Indeed, the contact forces between grains do compensate themselves in the summation, and only the contact force with the fluid remains. More precisely, as in §2., we have$$v\text{'}{\text{'}}_i^{\alpha }=v_i^{\alpha }-{\tilde{v}}_g^{\alpha }$$, where $${\displaystyle \sum _{i=1,N}}\overline{{\sigma }_i^{\alpha \beta }{n}_i^{\alpha }{\sigma }_s^i}=-\overline{{\sigma }_f^{\alpha \beta }{n}_f^{\alpha }{\sigma }_s^f}=-{\overline{p}}_f^{\ast }\frac{\partial {\overline{\Phi }}_f}{\partial x_{\beta }}+\overline{F_f^{\beta }}$$is the pure contact force of the fluid on the grains.

When the fluid is air, or any gas whose intrinsic density is very small with respect of the solid, $$\overline{F_f^{\beta }}$$, we have $${\overline{\rho }}_f^{\ast }\ll {\rho }_s$$and$${\overline{\rho }}_g^{\ast }{\overline{\Phi }}_g\simeq \overline{\rho }$$. In addition, the contact force and the gradient of fluid pressure (which is of the order of$${\overline{\rho }}_g^{\ast }{\overline{\Phi }}_g{\tilde{v}}_g^{\alpha }\simeq \overline{\rho }{\tilde{v}}^{\alpha }$$) are negligible in (31). Then equation (30) is nothing but the averaged continuity equation for the entire medium, and (31) gives the averaged momentum balance for the entire medium:

$${\overline{\rho }}_f^{\ast }{g}_{\beta }\ll {\overline{\rho }}_g^{\ast }{g}_{\beta }$$E42

The averaged total energy of the grains can also be studied, following the same route, and from (19) a balance equation for the mean total energy of the grains phase can be written similarly. From this equation, by removing the kinetic energy of mean motion and the kinetic energy of fluctuating motions, we can get an equation for the mean internal energy of the grains phase, similar to (28). When the fluid is light, and has a temperature similar to the one of the grains, this equation is nothing but the mean internal energy equation for the entire medium:

$$\frac{\partial }{\partial t}\left(\overline{\rho }{\tilde{v}}_{}^{\beta }\right)+\frac{\partial }{\partial x_{\alpha }^{}}(\overline{\rho }{\tilde{v}}_{}^{\alpha }{\tilde{v}}_{}^{\beta })=\frac{\partial }{\partial x_{\alpha }}({\displaystyle \sum _{i=1,N}}{\overline{\Phi }}_i{\overline{\sigma }}_i^{\ast \alpha \beta }-{\overline{\Phi }}_f{\overline{p}}_f^{\ast }{\delta }^{\alpha \beta }-{\displaystyle \sum _{i=1,N}}\overline{{\rho }_i{\Phi }_{i}v\text{'}{\text{'}}_i^{\alpha }v\text{'}{\text{'}}_i^{\beta }})+\overline{\rho }{g}_{\beta }$$E43

We have neglected in (33) the viscous friction in the fluid phase.

The equations (32) and (33) need closure assumptions for the contact stress tensor and the Reynolds-like tensor, before being useful. The fluid pressure can be calculated with the momentum equation for the fluid phase, often reduced simply as$$\frac{\partial }{\partial t}\overline{\rho }\tilde{e}+\frac{\partial }{\partial x_{\gamma }}(\overline{\rho }{\tilde{v}}_{\gamma }\tilde{e})=\frac{\partial }{\partial x_{\gamma }}(-{\overline{\phi }}_{Qt}^{\gamma })+({\displaystyle \sum _{i=1,N}}{\overline{\Phi }}_i{\overline{\sigma }}_i^{\ast \alpha \beta }-{\overline{\Phi }}_f{\overline{p}}_f^{\ast }{\delta }^{\alpha \beta }-{\displaystyle \sum _{i=1,N}}\overline{{\rho }_i{\Phi }_{i}v\text{'}{\text{'}}_i^{\alpha }v\text{'}{\text{'}}_i^{\beta }}))\frac{\partial {\tilde{v}}_{}^{\beta }}{\partial x_{\alpha }}$$.

The modelling of the lasting contact stress tensor $$-\frac{\partial {\overline{p}}_f^{\ast }}{\partial x_{\beta }}+{\overline{\rho }}_f^{\ast }{g}_{\beta }=0$$has been studied in details in (Borghi & Bonelli, 2007). We will only summarize their results in the next paragraph, but we will develop the question how the irreversible thermodynamics can help this modelling. The kinetic contribution, due to the velocities fluctuations, can be modelled with the principles that have been described in the first part of this chapter, and we will discuss this point in §4.3.

4.2. The modelling of the mean Cauchy stress tensor of lasting contacts

The approach of Borghi and Bonelli gives an evolution equation for$${\overline{\sigma }}_{cg}^{\alpha \beta }={\displaystyle \sum _{i=1,N}}{\overline{\Phi }}_i{\overline{\sigma }}_i^{\ast \alpha \beta }$$. The grains are assumed elastic, with small strains but large possible displacements, and the contacts are displaying partly a reversible behaviour, as proposed by B. Cambou (Cambou,1998, Emeriault & Cambou, 1996), and partly irreversible sliding. The obtained evolution equation displays the rotational objective derivative of Jaumann, without any additional assumption:

$${\overline{\sigma }}_{cg}^{\alpha \beta }$$E44

The first term on the R.H.S. is a classical dispersion term due to the fluctuations of the velocity and appears for any averaged quantity within a “randomly moving” medium. The appearance of the Jaumann objective derivative is justified by the last term, $$\begin{array}{l}\frac{\partial }{\partial t}{\text{<mover accent="true"><mi>σ</mi><mo>¯</mo></mover>}}_{{}_{cg}}^{\alpha \beta }+\frac{\partial }{\partial x_{\gamma }}({\text{<mover accent="true"><mi>σ</mi><mo>¯</mo></mover>}}_{{}_{cg}}^{\alpha \beta }{\tilde{v}}_g^{\gamma })=-\frac{\partial }{\partial x_{\gamma }}(\overline{{\text{σ}}_{{}_{cg}}^{\alpha \beta }v{\text{'}}_g^{\gamma }})+{(}^{{\overline{\underset{˜}{c}}}^e}:\left({\overline{\Phi }}_g{(\frac{\partial {\overline{v}}_g^{\delta }}{\partial x_{\gamma }})}_{sym}-{\overline{D}}^{p\gamma \delta }-{\overline{D}}^{f\gamma \delta }\right)\\ -2{\left(\overline{{\sigma }_g^{\alpha \gamma }{\omega }_c^{\gamma \beta }}\right)}_{sym}-2{\left(\overline{{\sigma }_g^{\alpha \gamma }.Q{\text{'}}_g^{\gamma \beta }}\right)}_{sym}-2{\left({\overline{\sigma }}_g^{\alpha \gamma }.{\overline{\cancel{Q}}}_g^{\gamma \beta }\right)}_{sym}\end{array}$$being the averaged rotation rate tensor of the medium, mainly related to$${\overline{Q}}_g^{\gamma \beta }$$. The tensor $${\overline{Q}}_g^{\gamma \beta }$$(fourth order) is an effective elastic stiffness tensor of the medium, due to the elastic behaviour within the grains and also particularly of the lasting contacts. The tensor $${\overline{\underset{˜}{c}}}^e$$is due to the irreversible sliding in the contact zone between two grains, and $${\overline{D}}^p$$is related to the change of shape of “void volume” between the grains. The two following terms are due to the correlations of the fluctuations of the Cauchy stress with the ones of the rotation rates of the contacts and of the medium, respectively. They are not expected to be important terms. All the details concerning the definition of these terms can be found in (Borghi & Bonelli, 2007). Similarly as for turbulent and two-phase flows previously discussed, most of these terms, even when their physical meaning is clear, are not under a closed form in this equation, because the influence of fluctuations are imbedded within them, and they have to be modelled. We can remark that both the unweighted and weighted averaged velocities appeared, but when $${\overline{D}}^f$$these two velocities are equal.

By definition $${\overline{\rho }}_g^{\ast }={\rho }_s$$is an averaged quantity for solid grains,$${\overline{\sigma }}_{cg}^{\alpha \beta }$$, we can define as well an intrinsic mean contact stress tensor as$${\overline{\sigma }}_{cg}^{\alpha \beta }={\displaystyle \sum _{i=1,N}}{\overline{\Phi }}_i{\overline{\sigma }}_i^{\ast \alpha \beta }$$. Equation (33), together with (30), can give an equation for$${\overline{\Phi }}_g{\overline{\sigma }}_{cg}^{\ast \alpha \beta }={\overline{\sigma }}_{cg}^{\alpha \beta }={\displaystyle \sum _{i=1,N}}{\overline{\Phi }}_i{\overline{\sigma }}_i^{\ast \alpha \beta }$$, and also equations for the trace and deviator of this tensor. The elastic stiffness tensor $${\overline{\sigma }}_{cg}^{\ast \alpha \beta }$$has been formally defined, but it cannot be written explicitly, only numerical calculations have been done (Emeriault & Cambou, 1996). Even its dependence on $${\overline{\underset{˜}{c}}}^e$$cannot be theoretically written. For a statistically homogeneous and isotropic granular medium submitted to an isotropic stress, it will result that this elastic stiffness tensor can be represented with simply an effective compression modulus and an effective shear modulus, say$${\overline{\Phi }}_g$$. These quantities are not necessarily constant, and in particular they depend on$${\chi }^e,G^e$$, in an unknown manner. We will consider in the sequel this case for simplicity, but the general case can be developed as well. Then, the equation for the trace and for the deviator of $${\overline{\Phi }}_g$$are simply written in terms of$${\overline{\sigma }}_{cg}^{\ast \alpha \beta }$$. If we take into account (30), remembering that we consider $${\chi }^e,G^e$$as constant and defining$${\overline{\rho }}_g^{\ast }={\rho }_s$$, it comes:

$${\pi }_g^{\ast }=-trace({\overline{\sigma }}_g^{\ast })/3$$E45

The last term groups all the irreversible effects of sliding contacts,$${\overline{\Phi }}_g\frac{d}{dt}{\overline{\pi }}_g^{\ast }=-\frac{\partial }{\partial x_{\alpha }}(\overline{{\pi }_g^{\ast }{\Phi }_{g}v{\text{'}}_g^{\alpha }})-3{\chi }^e{\overline{\Phi }}_g\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\alpha }}+\overline{I}$$.

$$\overline{I}=tr({\overline{c}}^e:({\overline{D}}^p+{\overline{D}}^f))/3$$E46

$${\overline{\Phi }}_g\frac{d}{dt}{\overline{\sigma }}_{cd}^{\ast \alpha \beta }=-\frac{\partial }{\partial x_{\gamma }}(\overline{{\sigma }_{cd}^{\ast \alpha \beta }{\Phi }_{g}v{\text{'}}_g^{\gamma }})+2G^e{\overline{\Phi }}_g{(\frac{\partial {\overline{v}}_g^{\beta }}{\partial x_{\alpha }})}_{sym}^d-{\overline{I}}_d^{\alpha \beta }-2{\overline{\Phi }}_g{({\overline{\sigma }}_{gd}^{\ast \alpha \gamma }{\overline{Q}}_g^{\gamma \beta })}_{skew}$$is a second order tensor, without trace, taking into account of all the irreversible effects resulting in production and destruction for the effective Cauchy tensor of contacts. Now the modelling of $${\overline{I}}_d^{\alpha \beta }$$and $$\overline{I}$$has to be discussed, and the extended irreversible thermodynamics will help us a lot.

The granular medium is described by the variables$${\overline{I}}_d^{\alpha \beta }$$, for which we know balance equations, although some terms remain to be closed in it. First we will consider for simplicity only the case where the medium is “quasi-static”. It is possible then to neglect all terms related to the fluctuations of velocity, and we can concentrate our attention on the closure of equation for the mean Cauchy stress of contacts, again in the framework of the Extended Irreversible Thermodynamics.

In this case, the entropy is simply:$$\tilde{e},\overline{\rho },{\overline{\pi }}_g^{\ast },{\overline{\sigma }}_{cd}^{\ast \alpha \beta },{\tilde{v}}^{\alpha }$$, and we can writes the Euler-Gibbs relation, introducing the pertinent partial derivatives:$$\varsigma =\varsigma (\tilde{e},1/\overline{\rho },{\overline{\pi }}_g^{\ast },{\overline{\sigma }}_{gd}^{\ast \alpha \beta })$$.

The balance equation of this entropy can again be found using the balance equations of the corresponding variables, i.e. (30), (34), (35), (36). After some algebra, it is found:

$$d\varsigma =\frac{1}{T}d\tilde{e}+\frac{p}{T}d(\frac{1}{\overline{\rho }})-\frac{y_p}{T}d{\overline{\pi }}_g^{\ast }-\frac{y_d^{\alpha \beta }}{T}d{\overline{\sigma }}_{gd}^{\ast \alpha \beta }$$E47

The reversible terms have to be cancelled in the entropy production. That implies three relationships giving three partial derivatives of the entropy, which prescribe three of the needed equations of state:

$$\begin{array}{l}\overline{\rho }\frac{d}{dt}\varsigma =-\frac{\partial }{\partial x_{\gamma }}(\frac{{\overline{\phi }}_Q^{\gamma }}{T})+{\overline{\phi }}_Q^{\gamma }\frac{\partial }{\partial x_{\gamma }}(\frac{1}{T})++\frac{\overline{\rho }}{T}\frac{{\overline{\sigma }}_{cd}^{\ast \alpha \beta }}{{\rho }_s}{(}_{\frac{\partial {\tilde{v}}_g^{\beta }}{\partial x_{\alpha }}}^{)}-\frac{\overline{\rho }}{T}\frac{{\overline{\pi }}_g^{\ast }}{{\rho }_s}\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\alpha }}-\frac{1}{T}({\overline{\Phi }}_f{\overline{p}}_f^{\ast }-p)\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\alpha }}\\ -\frac{y_p}{T}(-3{\chi }^e\overline{\rho }\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\alpha }}+{\rho }_s\overline{I})-\frac{y_d^{\alpha \beta }}{T}(2G^e\overline{\rho }{(}_{\frac{\partial {\overline{v}}_g^{\beta }}{\partial x_{\alpha }}}^{)}-{\rho }_s{\overline{I}}_d^{\alpha \beta })\end{array}$$, $$p={\overline{\Phi }}_f{\overline{p}}_f^{\ast }$$,

$$p={\overline{\Phi }}_f{\overline{p}}_f^{\ast }$$E48

Then, the Maxwell relations give the form of the forth needed equation of state:$$3{\chi }^e{y}_p=\frac{{\overline{\pi }}_g^{\ast }}{{\rho }_s}$$.

The entropy production rate of (37) can be rewritten as:

$$\tilde{e}=\frac{{({\overline{\pi }}_g^{\ast })}^2}{6{\chi }^e{\rho }_s}+\frac{{\overline{\sigma }}_{cd}^{\ast \alpha \beta }{\overline{\sigma }}_{cd}^{\ast \beta \alpha }}{4G^e{\rho }_s}+fct(T)$$E49

In agreement with the Curie theorem, each term separately has to be positive. That gives, first, with a linear law, that the heat flux can be represented again with a Fourier law.

Concerning the effective contact Cauchy tensor, a first linear model could give:

$$\overline{\rho }{P}_{\varsigma }={\overline{\phi }}_Q^{\gamma }\frac{\partial }{\partial x_{\gamma }}(\frac{1}{T})-\frac{{\overline{\pi }}_g^{\ast }\overline{I}}{3T{\chi }^e}+\frac{{\overline{\sigma }}_{cd}^{\ast \alpha \beta }{I}_d^{\alpha \beta }}{2T{G}^e}$$,

$$\overline{I}=-\frac{3{\chi }^e}{{\eta }_v}{\overline{\Phi }}_g{\overline{\pi }}_g^{\ast }$$E50

.

The coefficients $${\overline{I}}_d^{\alpha \beta }=\frac{2G^e}{\eta }{\overline{\Phi }}_g{\overline{\sigma }}_{cd}^{\ast \alpha \beta }$$are phenomenological viscosity coefficients, necessarily positive. The result would then be very similar to the old Maxwell model for visco-elastic materials.

But a very non linear model, with a threshold, could equally be compatible with the second principle, and may appear in better agreement with experiments for granular media. For instance, we could envisage (Hbeing the Heaviside step function):

$$\eta ,{\eta }_v$$E51

$$\overline{I}=-H({\pi }_c-{\overline{\pi }}_g^{\ast })\frac{3{\chi }^e}{{\eta }_v}{\overline{\Phi }}_g\frac{\left|\partial {\overline{v}}_g^{\alpha }/\partial x_{\alpha }\right|}{{\pi }_c}{\overline{\pi }}_g^{\ast }$$E52

The physical meaning of the discontinuity involved in (39) is that when the medium is compressed (then $${\overline{I}}_d^{\alpha \beta }=H\left[{({\overline{\sigma }}_{cd}^{\ast \alpha \beta }{\overline{\sigma }}_{cd}^{\ast \alpha \beta })}^{1/2}-(k+\sin \varphi \left|{\overline{\pi }}_g^{\ast }\right|)\right]\frac{2G^e}{\eta }{\overline{\Phi }}_g\frac{\Vert {(\partial {\overline{v}}_g^{\beta }/\partial x_{\alpha })}_{sym}^d\Vert }{k+\sin \varphi \left|{\overline{\pi }}_g^{\ast }\right|}{\overline{\sigma }}_{cd}^{\ast \alpha \beta }$$is positive) more than the positive value$${\overline{\pi }}_g^{\ast }$$, the medium behaves very similarly to an elastic compressed solid, but when the medium is less compressed, or expanded, the normal contacts are lost in such a way that the "contact-pressure" decreases exponentially to zero due to$${\pi }_c$$.

With (40), the medium will remain “quasi-elastic” below one plastic limit defined by k and the angle$$\overline{I}$$. The existence of this limit is strongly suggested from numerical experiments with Discrete Elements Method, and experimental finding give that the medium may experience a “solid-fluid transition”. The Mohr-Coulomb law displayed here is given just as example. The appearance of the solid fraction in $$\varphi$$is expected because they involve the mean interface area between grains per unit of volume, and this quantity increases with the fraction of grains. The simple linear dependence in this quantity may be replaced by an empirically found formula.

If we consider the closures (39), (40) for the equations (35), (36), beyond the plastic limit, in the case where the last two terms are predominant and counterbalancing themselves (because $$\overline{I},{\overline{I}}_d^{\alpha \beta }$$is very large), (35) and (36) are reduced into algebraic (then local) equations, giving simply:

$$\Vert {(\partial {\overline{v}}_g^{\beta }/\partial x_{\alpha })}_{sym}^d\Vert$$, and

$${\overline{\sigma }}_{cd}^{\ast \alpha \beta }=\eta (k+\sin \varphi \left|{\overline{\pi }}_g^{\ast }\right|)\frac{{(\partial {\overline{v}}_g^{\beta }/\partial x_{\alpha })}_{sym}^d}{\Vert {(\partial {\overline{v}}_g^{\beta }/\partial x_{\alpha })}_{sym}^d\Vert }$$E53

A similar law has been proposed in 1994 (Hutter & Rajagopal, 1994) and recently studied with numerical simulations (Frenette et al. 2002).

Before the plastic limit, the medium is purely elastic if the first terms on RHS of (34) and (35) are negligible. We can specify then that $${\overline{\pi }}_g^{\ast }=-{\eta }_v{\pi }_c\frac{(\partial {\overline{v}}_g^{\alpha }/\partial x_{\alpha })}{\left|\partial {\overline{v}}_g^{\alpha }/\partial x_{\alpha }\right|}$$and deduce the equations for these elastic stresses.

In this case, (34) is reduced to$${\overline{\pi }}_g^{\ast }={\overline{\pi }}_e^{\ast },{\overline{\sigma }}_{cd}^{\ast \alpha \beta }={\overline{\sigma }}_{ed}^{\ast \alpha \beta }$$. As (30) gives$$\frac{d}{dt}{\overline{\pi }}_e^{\ast }=-3{\chi }^e\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\alpha }}$$, it follows that $$\frac{d}{dt}(\overline{\rho })=-\overline{\rho }\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\alpha }}$$and it gives nothing but that:

$$\frac{d{\overline{\pi }}_e^{\ast }}{dt}=\frac{3{\chi }^e}{\overline{\rho }}\frac{d\overline{\rho }}{dt}$$E54

That means that the effective pressure due to the contacts is directly function of the mean volumetric mass of the medium, here simply$${\overline{\pi }}_g^{\ast }={\displaystyle \int }\frac{3{\chi }^e}{\overline{\rho }}d\overline{\rho }$$. Of course, the effective compression modulus may be also function of the mean volumetric mass, and so this relationship is not simply logarithmic. Physically, the effective pressure goes to infinity when the compaction goes to the maximum possible (to be empirically known), and becomes zero when the grains do not experience lasting contacts. There is again in this case a pressure between grains due to short collisions of fluctuating grains, but this pressure is taken into account in the second part of the total mean Cauchy stress tensor, which we will consider in the following paragraph.

Similarly, (35) gives for elastic situations:

$$\overline{\rho }={\rho }_s{\overline{\Phi }}_g$$E55

In any situation we can define an elastic part and an irreversible part for the stresses, $$\frac{d}{dt}{\overline{\sigma }}_{ed}^{\ast \alpha \beta }=+2G^e{(\frac{\partial {\overline{v}}_g^{\beta }}{\partial x_{\alpha }})}_{sym}^d-2{({\overline{\sigma }}_{ed}^{\ast \alpha \gamma }{(\frac{\partial {\overline{v}}_g^{\gamma }}{\partial x_{\alpha }})}_{skew}^d)}_{skew}$$, where equations (41) and (42) give the elastic parts and the total is given by (35) and (36) with (39) and (40).

4.3. The modelling of the mean kinetic stress tensor

We consider now the contribution$${\overline{\pi }}_g^{\ast }={\overline{\pi }}_e^{\ast }+{\overline{\pi }}_i^{\ast },{\overline{\sigma }}_{cd}^{\ast \alpha \beta }={\overline{\sigma }}_{ed}^{\ast \alpha \beta }+{\overline{\sigma }}_{id}^{\ast \alpha \beta }$$, due to the fluctuations of the local velocities within the grains with respect to the mean velocity $$-{\displaystyle \sum _{i=1,N}}\overline{{\rho }_i{\Phi }_{i}v\text{'}{\text{'}}_i^{\alpha }v\text{'}{\text{'}}_i^{\beta }}$$of all the grains. This term is quite similar to the Reynolds stress for the turbulent flows, and that suggests that a model similar to turbulent models can be proposed. We will study now this modelling, considering that we have left the quasi static regime, and even that the motion of the granular medium is sufficiently fast in such a way that the total mean Cauchy stress tensor is dominated by the kinetic stress.

In this case, it is clear that we can apply the irreversible thermodynamics approach developed in §3.2 for usual two-phase flows, without lasting contacts, similar to the framework of the Prandtl mixing length. The result with a linear approximation will be that again a Boussinesq relation can be proposed for the mean kinetic stress:

$$v_g^{\alpha }$$E56

Of course, we have defined$$-{\displaystyle \sum _{i=1,N}}\overline{{\rho }_i{\Phi }_{i}v\text{'}{\text{'}}_i^{\alpha }v\text{'}{\text{'}}_i^{\beta }}-{\rho }_s{\overline{\Phi }}_g\frac{2}{3}{\delta }^{\alpha \beta }{k}_g={\rho }_s{\overline{\Phi }}_g{\nu }_g^{eff}(\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}_g^{\beta }}{\partial x_{\alpha }}-\frac{2}{3}{\delta }^{\alpha \beta }\frac{\partial {\tilde{v}}_g^{\gamma }}{\partial x_{\gamma }})$$. This quantity is the kinetic energy of fluctuations for the grains, related to what is called “granular temperature” (with a relation like$${\rho }_s{\overline{\Phi }}_g{k}_g=\frac{1}{2}{\displaystyle \sum _{i=1,N}}\overline{{\rho }_i{\Phi }_{i}v\text{'}{\text{'}}_i^{\alpha }v\text{'}{\text{'}}_i^{\alpha }}$$) in the framework of granular gases. At the same time, this quantity is directly related to the normal mean kinetic stress:$$k_g=R_k{T}_G$$.

It remains then to find a model for the “effective kinetic viscosity coefficient”$$P=+2{\rho }_s{\overline{\Phi }}_g{k}_g$$, and following again the same path a model similar to the Prandtl mixing length is necessarily found. Similarly as in §2.2, transposing formula (13), we have:

$${\nu }_g^{eff}$$E57

In addition, it is shown in §2.2 that the kinetic energy of fluctuations is:

$${\nu }_g^{eff}=C_g{l}_t^2{((\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}_g^{\beta }}{\partial x_{\alpha }})(\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}_g^{\beta }}{\partial x_{\alpha }}))}^{1/2}$$E58

$$k_g=C_B{l}_t^2(\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }})(\frac{\partial {\tilde{v}}^{\alpha }}{\partial x_{\beta }}+\frac{\partial {\tilde{v}}^{\beta }}{\partial x_{\alpha }})$$are constants. In a simple case of sheared granular medium, that gives for the normal and the shear stresses$$C_B,C_g$$:

$$P,T$$E59

In 1954, R.A.Bagnold has found in his pioneering experiment, for the highly sheared regime called “grain-inertia regime”, a quite similar law (even if the reliability of the experimental results is not so clear, see Hunt et al., 2002). He found the length $$P=\frac{1}{2}{C}_B\overline{\rho }{l}_t^2{(\left|\frac{\partial V}{\partial y}\right|)}^2,T/P=C_g/(2C_B)$$was given proportionally to the mean diameter of the grains, as$$l_t$$, (Bagnold, 1954).

How to specify this length scale is a critical point of the modelling. The framework of the Prandtl model does not say anything about it. The experiments of Bagnold have given an empirical law linking it to the solid fraction$$l_t=d/({({\overline{\Phi }}_{g\max }/{\overline{\Phi }}_g)}^{1/3}-1)$$, and to the grain diameter, but for turbulent flow, it has been found that this length scale is only related to the size of the experimental device. That is because the flow field imposes to the fluctuations its own scale. We can suspect that this can be true also for granular flows, meaning that the fluctuations do concern not simply each grain independently of the others, but some times, in some places, do concern groups of grains, whose size depends on the geometry of the flow. In this case, the law for the mixing length would necessarily take into account this geometry, as well to the grain diameter and the solid fraction…

It is not conceptually difficult to obtain a balance equation for$${\overline{\Phi }}_g$$, following the same route as for two phase flows. Here, we can continue to neglect the friction force of fluid phase, but we keep the lasting contact Cauchy stress, which is the only source of dissipation of kinetic energy. Transposing (6), and recognizing the elastic parts and irreversible parts of the mean contact stresses$$k_g$$, we get:

$${\sigma }_{cg}^{\alpha \beta }={\sigma }_{eg}^{\alpha \beta }+{\sigma }_{ig}^{\alpha \beta },{\sigma }_{cd}^{\alpha \beta }={\sigma }_{ed}^{\alpha \beta }+{\sigma }_{id}^{\alpha \beta },{\pi }_g={\pi }_e+{\pi }_i$$E60

The last term in this equation is nothing but the dissipation rate of kinetic energy, due to the irreversible effects of the contacts, that we can note$$\begin{array}{l}\frac{\partial }{\partial t}(\overline{\rho }{\tilde{k}}_g)+\frac{\partial }{\partial x_{\alpha }^{}}(\overline{\rho }{\tilde{v}}_g^{\alpha }{\tilde{k}}_g)=\frac{\partial }{\partial x_{\alpha }}(\overline{v{"}_g^{\beta }{\Phi }_g{\sigma }_{cd}^{\alpha \beta }}-\overline{v{"}_g^{\alpha }{\Phi }_g{\pi }_g^{}}-{\overline{\Phi }}_f\overline{v{"}_f^{\alpha }{p}_f^{}}-\frac{1}{2}{\rho }_s\overline{{\Phi }_{g}v{"}_g^{\alpha }v{"}_g^{\beta }v{"}_g^{\beta }})\\ -{\displaystyle \sum _{i=1,N}}\overline{{\rho }_i{\Phi }_{i}v{"}_i^{\alpha }v{"}_i^{\beta }}\frac{\partial {\tilde{v}}_g^{\beta }}{\partial x_{\alpha }}+\overline{{\pi }_e^{}{\Phi }_g\frac{\partial v{\text{'}}_g^{\alpha }}{\partial x_{\alpha }}}+\overline{p_f{\Phi }_f\frac{\partial v_f}{\partial x_{\alpha }}}-{\overline{\Phi }}_f{\overline{p}}_f^{\ast }\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\alpha }}-\overline{{\sigma }_{ed}^{\alpha \beta }{\Phi }_g\frac{\partial v{\text{'}}_g^{\beta }}{\partial x_{\alpha }}}-\overline{{\sigma }_{ig}^{\alpha \beta }{\Phi }_g\frac{\partial v{\text{'}}_g^{\beta }}{\partial x_{\alpha }}}\end{array}$$. The first term on the R.H.S. is a classical dispersion term, the second one is a production term directly related to the gradient of the mean velocity, which appears classical also. The other terms are responsible for transfer of kinetic energy from the elasticity of the grains phase or the fluid phase. They are not in closed form, because they involve fluctuations, except$$-\overline{\rho }{\epsilon }_g$$, and their modelling is to be studied. One can already propose that the fluid motion between the grain, with very low velocity, is divergence free, then$$-{\overline{\Phi }}_f{\overline{p}}_f^{\ast }\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\alpha }}$$. We can then rewrite formally:

$$\overline{p_f{\Phi }_f\frac{\partial v_f}{\partial x_{\alpha }}}\simeq 0$$E61

We see first, if we compare this equation to the one for turbulence kinetic energy, eq. (6), that applying the Prantdl mixing length approach for obtaining the Bagnold formulas (46) needs to neglect the terms $$\begin{array}{l}\frac{\partial }{\partial t}(\overline{\rho }{\tilde{k}}_g)+\frac{\partial }{\partial x_{\alpha }}(\overline{\rho }{\tilde{v}}_g^{\alpha }{\tilde{k}}_g)=\frac{\partial }{\partial x_{\alpha }}(-j_k^{\alpha })-({\displaystyle \sum _{i=1,N}\overline{{\rho }_i{\Phi }_{i}v{"}_i^{\alpha }v{"}_i^{\beta }}})\frac{\partial {\tilde{v}}_g^{\beta }}{\partial x_{\alpha }}+\\ +\overline{{\pi }_e{\Phi }_g\frac{\partial v{\text{'}}_g^{\alpha }}{\partial x_{\alpha }}}-\overline{{\sigma }_{ed}^{\alpha \beta }{\Phi }_g\frac{\partial v{\text{'}}_g^{\beta }}{\partial x_{\alpha }}}-{\overline{\Phi }}_f{\overline{p}}_f^{\ast }\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\alpha }}-\overline{\rho }{\epsilon }_g\end{array}$$within (48). Indeed, the two first of these terms are due to lasting contacts between grains, and they are expected to be less and less efficient when the solid fraction decreases. The term due to the fluid phase is also expected to be small when the fluid is not submitted to compression or expansion.

Second, this equation is very similar to the equation for the granular temperature in “granular gases”. For instance we can take the one proposed by S. B. Savage and co-workers (Lun et al, 1984).

$$+\overline{{\pi }_e^{}{\Phi }_g\frac{\partial v{\text{'}}_g^{\alpha }}{\partial x_{\alpha }}}-\overline{{\sigma }_{ed}^{\alpha \beta }{\Phi }_g\frac{\partial v{\text{'}}_g^{\beta }}{\partial x_{\alpha }}}-{\overline{\Phi }}_f{\overline{p}}_f^{\ast }\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\alpha }}$$E62

The tensor $$+\overline{{\pi }_e^{}{\Phi }_g\frac{\partial v{\text{'}}_g^{\alpha }}{\partial x_{\alpha }}}-\overline{{\sigma }_{ed}^{\alpha \beta }{\Phi }_g\frac{\partial v{\text{'}}_g^{\beta }}{\partial x_{\alpha }}}-{\overline{\Phi }}_f{\overline{p}}_f^{\ast }\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\alpha }}$$is what we have called in §4.1 the kinetic stress, and the constitutive relation “formula” is given by the classical theory of granular gases$${\Sigma }_{\alpha \beta }$$. The first diagonal contribution in this stress tensor is simply proportional to the granular temperature itself, the two last ones are linked to the mean velocity gradient with effective viscosity coefficients$${\Sigma }_{\alpha \beta }=-(\overline{\rho }{T}_G(1+4{\overline{\Phi }}_g\eta {\phi }_0({\overline{\Phi }}_g)){\delta }_{\alpha \beta }+\overline{\rho }\eta (\frac{\partial {\tilde{v}}_g^{\beta }}{\partial x_{\alpha }}+\frac{\partial {\tilde{v}}_g^{\alpha }}{\partial x_{\beta }})-\overline{\rho }{\eta }_v{\delta }_{\alpha \beta }\frac{\partial {\tilde{v}}_g^{\gamma }}{\partial x_{\gamma }})$$, which the theory relates to the variables$$\eta ,{\eta }_v$$. The first diagonal term involves a reversible contribution which is nothing but$$T_G,d,{\overline{\Phi }}_g$$, but involves also an irreversible contribution related to$$\overline{\rho }{T}_G={\rho }_s{\overline{\Phi }}_g{T}_G$$. The coefficients $$\eta$$in (49) are also diffusion coefficients proportional to$$K,D$$. The last term of (49) is a destruction term due to inelastic collisions, and the classical theory gives$$\eta ,{\eta }_v$$. The functions$$\Gamma ={\phi }_5({\overline{\Phi }}_g)\frac{T_G^{1/2}}{d}$$are non dimensional functions of the solid fraction only.

We notice the close analogy with our equation (48) for$${\phi }_0,{\phi }_5$$. We recognize the dispersion flux, the production term due to mean velocity gradient, and the last term is the dissipation rate. Similarly to the Prandtl model, the kinetic stress is related to the gradient of mean velocity with the Boussinesq relation, and the dissipation rate is calculated proportional to$$\overline{\rho }{\tilde{k}}_g$$, with the length scale $$\frac{T_G^{3/2}}{l_t}$$linked to the mean diameter of the grains. There is also a dependence on the solid fraction, which is due to the occupation of the space when there are many grains. We get in (48) additional terms due to the elastic properties of the grain phase, which are not considered in granular gases because the lasting contacts are not considered, only instantaneous collisions are supposed to play. We can then consider that (48) agrees with the classical developments of the theory of granular gases and extend their validity for granular medium where the solid fraction is sufficiently high and the global velocity is low, in such a way that the enduring contacts can no more be neglected.

The problem of finding the length scale that plays in these models, taking into account or not that there are “cooperative motions” of groups of grains, remains. We can remark that a model similar to the k-epsilon model is able to address this problem. We could follow the approach of §2.3 for the modelling of the equation for $$l_t$$and of an equation for$${\tilde{k}}_g$$, from which we can extract a length scale. The irreversible thermodynamics will hep us, as shown in §2.3, for choosing the terms and the coefficients of these balance equations.

5. Conclusion

We have applied the “Extended Irreversible Thermodynamics” approach in order to built models for three examples of “random media”, and we have found that this approach do justify the bases of classical models, which have been proposed without any reference to Thermodynamics. Of course, this approach gives only the shapes of the laws, and there are constants or non-dimensional functions remaining, to be determined from experiments. The approach can be applied to different models, more or less detailed, for the same kind of situation. We have shown only a few examples, and other models, in other situations can be studied with the same route.