# Modelling of relative velocity, velocity fluctuations and their interactions for two-fluid models by Stationary Action Principle

Ward Haegeman, Giuseppe Orlando, Samuel Kokh, Marc Massot

► To cite this version:

Ward Haegeman, Giuseppe Orlando, Samuel Kokh, Marc Massot. Modelling of relative velocity, velocity fluctuations and their interactions for two-fluid models by Stationary Action Principle. 2025. [hal-05224712](https://hal.science/hal-05224712v1)

HAL Id: hal-05224712

https://hal.science/hal-05224712v1

Preprint submitted on 26 Aug 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
MODELLING OF RELATIVE VELOCITY, VELOCITY FLUCTUATIONS AND
THEIR INTERACTIONS FOR TWO-FLUID MODELS BY STATIONARY
ACTION PRINCIPLE *

WARD HÆGEMAN<fnref n="1" />,<fnref n="2" />, GIUSEPPE ORLANDO<fnref n="2" />, SAMUEL KOKH<fnref n="3" /> AND MARC MASSOT<fnref n="2" />

**Abstract.** The objective of this contribution is the derivation of a two-fluid model including a relative velocity between the two phases and velocity fluctuations, describing pseudo-turbulent effects, as internal variables based on Stationary Action Principle. The variational derivation, used to obtain the model, relies on the variation of a single trajectory related to the mass-weighted average velocity under the barotropic assumption. The model is hyperbolic, satisfies a second principle of thermodynamics, and admits either linearly degenerate or genuinely nonlinear characteristic fields. Moreover, the variational approach yields a fully closed model and its non-conservative products are uniquely defined for weak solutions in 1D, i.e. jump conditions can be derived. In the laminar case, when velocity fluctuations are negligible, we recover previously derived multi-fluid models which have been analyzed in several contributions. As such, the present framework allows for an original extension of the existing models to include velocity fluctuations of each phase for pseudo-turbulent flows, their coupling with the relative velocity between phases, as well as dissipative effects compatible with the thermodynamics of irreversible processes. Eventually, we provide a discussion of the limitations of the proposed model, especially regarding the extension to the open problem of non-barotropic flows.

**Résumé.** L'objectif de cette contribution est la dérivation d'un modèle bi-fluide incluant une vitesse relative entre les deux phases et les fluctuations de vitesse, décrivant des effets pseudo-turbulents, sous la forme de variables internes. Cette dérivation est basée sur le Principe de l'Action Stationnaire. La dérivation variationnelle du modèle repose sur la variation d'une trajectoire unique liée à la vitesse du centre de masse du mélange sous l'hypothèse barotropique. Le modèle est hyperbolique, satisfait un second principe de la thermodynamique et admet des champs caractéristiques linéairement dégénérés ou vraiment non linéaires. L'approche variationnelle fournit un modèle totalement fermé dont les produits non-conservatifs sont définis de manière non-ambiguë pour les solutions faibles dans le cas mono-dimensionnel dans le sens où des conditions de saut peuvent être dérivés. Dans le cas laminaire, lorsque les fluctuations de vitesse sont négligéables, nous retrouvons des modèles multi-fluides existants et analysés dans la littérature. Ainsi, le cadre variationnel permet une extension originale de ces modèles pour inclure des fluctuations de vitesse pour chaque phase dans le cas d'écoulements pseudo-turbulents et leur couplage avec la vitesse relative, ainsi que des effets dissipatifs compatibles avec la thermodynamique des processus irréversibles. Enfin, nous discutons des limites du modèle, en particulier en ce qui concerne l'extension au problème ouvert des écoulements non barotropes.

<footnote>* This work has been supported by the Agence Innovation Défense (AID), the CIEDS project OPEN-NUM-DEF (PI. L. Gouarin, T. Pichard, M. Massot), a PhD Grant from ONERA (W. Haegeman), as well as the HPC@Maths Initiative (PI. L. Gouarin and M. Massot), which we would like to acknowledge.</footnote>

<footnote n="1">1 DMPE, ONERA Université Paris Saclay, F-91123 Palaiseau, France</footnote>

<footnote n="2">2 CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris, Route de Saclay, Palaiseau, 91120, France</footnote>

<footnote n="3">3 Maison de la Simulation, UAR 3441, Bâtiment 565 PC 190 - Digiteo, CEA Saclay 91191, Gif-sur-Yvette, France</footnote>
# INTRODUCTION

The study of two-phase flows is of interest to many applications. Far from critical conditions, such flows may be modelled by two sets of Navier-Stokes equations on unsteady domains, coupled through appropriate jump conditions across the interface, which can be followed in time. Such descriptions are referred to as one-fluid models. Several numerical discretizations of these models have been developed, such as level-set methods, volume of fluid methods, and interface-tracking techniques. We refer, e.g., to [21] for an overview of interface-capturing methods for two-phase flows. However, for many applications, the multi-scale nature of the flow and the important topology changes that occur prohibit the use of such numerical strategies and an alternative description is needed. By allowing a description which does not require to explicitly reconstruct the interface, multi-fluid models are a popular alternative. The development of suitable models for compressible two-phase flows is still an active area of research, as no universally accepted model has yet been established, even when velocity fluctuations are neglected. Among the most widely used models for compressible two-phase flows are the Baer-Nunziato-type models [1], which assume that each phase evolves with its own pressure, temperature, and velocity, together with an evolution equation for the volume fraction of one of the two phases. Originally, the Baer-Nunziato (BN) model was proposed for detonation waves in granular explosives and then extended to liquid-gas two-phase flows in [29]. In the context of two-phase flows, several variants of the model exist as it relies on the introduction of some interfacial quantities which are the interfacial velocity and the interfacial pressure. Various closures for these quantities have been proposed in the literature, as they impact both the physics described and the mathematical properties of the resulting models. Moreover, the models are non-conservative, thus yielding some difficulties in the definition of weak solutions and jump conditions. Nevertheless, specific closures of the interface variables that guarantee the uniqueness of non-conservative products have been derived in [7].

Alternatively, compressible two-phase models based on the theory of thermodynamically compatible systems [6, 10, 11, 28] have been derived [24–26]. Similarly to the Baer-Nunziato models, they also allow for a complete description of non-equilibrium effects. In the barotropic setting, the models have been analyzed both from a theoretical and numerical point of view in the one-dimensional setting in [32, 33]. The formulation presented in [24–26] adopts an evolution equation for the mass-weighted average velocity and an evolution equation for the relative velocity between the two phases. Interestingly, in 1D, when recasted in a BN-type form, i.e. considering the evolution equations for the phasic momenta, the system turns out to be identical to the barotropic BN system with specific closures of the interfacial velocity and pressure [33]. In the multi-dimensional setting, lift forces related to the vorticity of the phasic velocities appear and modify the mathematical structure of the model.

The aim of this paper is to present a barotropic compressible two-fluid model that includes relative velocity and velocity fluctuations. We propose an original variational derivation based on the Stationary Action Principle relying on a single trajectory related to the mass-weighted average velocity. This approach allows to avoid the introduction of two family of trajectories as in [3, 9] and the difficulties associated to their interactions [8]. The derivation follows the generic Stationary Action Principle framework presented in [13]. Within this setting, the relative velocity between the two phases is an internal variable that contributes to the description of the state of the two-fluid mixture. The idea of considering the relative velocity as an internal variable has been already present in [24–26] and, previously, in [18, 19]. We extend this approach to enrich the set of state variables of the two-fluid mixture by including the velocity fluctuations as additional internal variables in order to describe pseudo-turbulent effects. This provides a novel way to include such effects in multi-fluid models.

In the laminar case, i.e. when velocity fluctuations are negligible, we recover the model proposed in [24–26]. However, the Stationary Action Principle framework provides a natural way to incorporate velocity fluctuations and to provide their coupling with the relative motion between the phases. The variational framework allows for the derivation of the reversible part of the flow. Dissipative source terms are then added compatibly with the second principle of thermodynamics. In the present work, we analyze the dissipation induced by drag which drives the system towards kinematic equilibrium. A reduced order model is derived through asymptotic
expansion near equilibrium, in which the relative velocity is given by a closed expression of the remaining flow variables, including the pseudo-turbulent energies.

The paper is structured as follows. In Section 1, we present the derivation of the model through Stationary Action Principle. We exhibit the link with the models obtained with the theory of thermodynamically compatible systems and provide a discussion about the issues and limitations for the extension to non-barotropic flows. In Section 2, the model derived in Section 1 is extended accounting for velocity fluctuations. Section 3 focuses on the inclusion of dissipative effects to ensure consistency with the second principle of thermodynamics and the derivation of reduced order models near equilibrium. Some conclusions and perspectives for the future developments conclude the manuscript.

# 1. TWO-FLUID MODEL WITH RELATIVE VELOCITY AS INTERNAL VARIABLE

In this section, we propose to use the Stationary Action Principle to derive the two-fluid model with relative velocity between the two phases, as velocity fluctuations will be accounted for in the following section. We first recall the Stationary Action Principle framework presented in [13].

## 1.1. The Stationary Action Principle framework

A generic framework to derive systems of conservation laws based on the Stationary Action Principle has been presented in [13]. We provide the main results of the framework as it will be used in Section 1.2 to derive the two-fluid model. Consider a Lagrangian density $\mathcal{L}$ that depends on the fluid velocity $\mathbf{u}$ and its density $\rho$. There are two sets of additional variables on which the Lagrangian depends, a set of transported variables $\mathcal{B} = \{b_i, i = 1..., n\}$ which satisfy $D_t b_i = 0$ — $D_t \bullet$ denoting the material derivative — and a set of free variables $\Xi = \{\xi_j, j = 1..., m\}$ and their derivatives. In particular, derivatives in time are expressed through material derivatives, as this choice is indeed more suitable to obtain a frame invariant expression of the resulting system. The general form of the Lagrangian density that we consider therefore reads

$$ \mathcal{L} \equiv \mathcal{L}(\mathbf{u}, \rho, (b_i)_{i=1}, (D_t \xi_j, \nabla \xi_j, \xi_j)_j). \quad (1) $$

The associated action $\mathcal{A}$ then corresponds to the integral of the Lagrangian density (1) over the entire space-time domain, i.e.

$$ \mathcal{A} = \int_t \int_{\mathcal{V}(t)} \mathcal{L} d\mathbf{x} dt. \quad (2) $$

The variational principle then dictates that the equations of motion are critical points of the Action, in the sense of Gateaux derivatives. However, in order to obtain the correct set of equations, an underlying Lagrangian description based on fluid particle trajectories is necessary and the Gateaux derivative is to be applied to these trajectories. For more details, we refer to [13] and references therein and without further ado we provide the generic system of equations derived from (1) which writes

$$ \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, \quad (3a) $$

$$ \partial_t (\rho b_i) + \nabla \cdot (\rho b_i \mathbf{u}) = 0, \quad i = 1, \dots, n, \quad (3b) $$

$$ \partial_t \mathbf{K} + \nabla \cdot (\mathbf{K} \otimes \mathbf{u}) + \nabla \cdot \mathbf{\Pi} = \mathbf{0}, \quad (3c) $$

$$ \partial_t \mathcal{M}_j + \nabla \cdot (\mathcal{M}_j \mathbf{u}) - \frac{\partial \mathcal{L}}{\partial \xi_j} + \nabla \cdot \mathbf{D}_j = 0, \quad j = 1, \dots, m. \quad (3d) $$

Here, $\mathbf{K}$ is the conjugate variable of the velocity and represents the generalized momentum. It obeys a conservation law in which $\mathbf{\Pi}$ is the generalized pressure tensor. The supplementary equations on the variables $\mathcal{M}_j$ result from the dependence of the Lagrangian to the material derivatives $D_t \xi_j$. The notations are defined as
follows:

$$ K^T = \frac{\partial \mathcal{L}}{\partial \mathbf{u}}, \quad \Pi = -\mathcal{L}^* \mathbf{Id} - \sum_{j=1}^{m} \nabla \xi_j \otimes D_j, \quad \mathcal{M}_j = \frac{\partial \mathcal{L}}{\partial(D_t \xi_j)}, \quad D_j^T = \frac{\partial \mathcal{L}}{\partial(\nabla \xi_j)}, \quad (4) $$

where $\mathbf{Id}$ is the identity tensor and $\mathcal{L}^* = \rho \frac{\partial \mathcal{L}}{\partial \rho} - \mathcal{L}$. The conventions used here are that vector fields are represented as column vectors. The derivative of a scalar field with respect to a vector being its differential, it is a row vector, its transpose is the scalar's gradient. Similarly, the gradient of a vector field is the transpose of its Jacobian matrix. Finally, the divergence of a matrix is the column vector obtained by taking the divergence of each of its rows.

One of the advantages of the variational approach is that system (3) admits the following supplementary conservation law

$$ \partial_t \mathcal{H} + \nabla \cdot \left( \mathcal{H} \mathbf{u} + \mathbf{\Pi}^T \mathbf{u} + \sum_{j=1}^{m} D_j D_t \xi_j \right) = 0, \quad (5) $$

where the Hamiltonian $\mathcal{H}$ is the Legendre transform with respect to $\mathbf{u}$ and the material derivatives of the free variables. As such, the Hamiltonian is a function of the prognostic variables $\mathbf{K}$ and $\mathcal{M}_j$ instead of $\mathbf{u}$ and $D_t \xi_j$. Its expression is given by

$$ \mathcal{H} = \mathbf{K}^T \mathbf{u} + \sum_{j=1}^{m} \mathcal{M}_j D_t \xi_j - \mathcal{L}. \quad (6) $$

Concerning the convexity requirement of the Legendre transform, in most physical systems, the Lagrangian is quadratic with respect to the velocity, and thus satisfies the convexity requirement. However, the variables $\xi_j$ may represent Lagrange multipliers used to enforce internal constraints. In this case, the Lagrangian may be linear in $D_t \xi_j$ and the convexity requirement for the Legendre transform therefore is not satisfied. However, their contribution cancels out as the Hamiltonian corresponds to the total energy and should only depend on the variables which bear a physical meaning.

## 1.2. Modelling relative velocity in two-fluid models

In order to derive a multi-fluid two-phase flow model which includes the relative velocity and its interactions with velocity fluctuations, we start with the derivation of a model including only relative velocity. Velocity fluctuations will be included afterwards in Section 2. One of the questions raised by multi-velocity models is that of interfacial velocity, i.e. what is the correct velocity $\mathbf{u}_I$ for the transport of the volume fraction? In the original Baer-Nunziato multi-fluid model [1], the volume fraction is transported at one of the phasic velocities, i.e. $\mathbf{u}_I = \mathbf{u}_k$, $k \in \{1, 2\}$. While this simple assumption may be acceptable in many applications, it may not suffice in flows with strong topology changes, for instance when the carrier phase in one region corresponds to the disperse phase in another region. A practical choice when aiming for an all-topology model is to consider the mass-weighted average velocity $\mathbf{u}_I = Y_1 \mathbf{u}_1 + Y_2 \mathbf{u}_2$, where $Y_k$, $k = 1, 2$ are the mass fractions which sum up to 1. This is the choice that we consider in this work. This velocity will be referred to as barycentric velocity in the following and, for the sake of convenience, we drop the index $I$, i.e. we set $\mathbf{u} = Y_1 \mathbf{u}_1 + Y_2 \mathbf{u}_2$.

The Stationary Action Principle being closely related to the Lagrangian representation, a two-velocity model would call for two families of trajectories, one for each phase. The interaction between the two corresponding reference frames in the variational principle is non-trivial [8]. Here, to avoid this difficulty, we choose to derive a model with a single reference frame attached to the barycentric velocity $\mathbf{u}$ and thus consider a single family of virtual trajectories<fnref n="1" /> associated to this velocity field. Velocity non-equilibrium is then described in this single

<footnote n="1">¹Here, the term virtual trajectory indicates that the trajectory, associated to the barycentric velocity, which is used in the variational principle, differs from the physical trajectories of the fluid particles.</footnote>
reference frame by an internal variable, which corresponds to the relative velocity $\mathbf{W} = \mathbf{u}_1 - \mathbf{u}_2$. Doing so, we write the total kinetic energy in terms of $\mathbf{u}$ and $\mathbf{W}$, i.e.

$$ \frac{1}{2}\alpha_1\rho_1|\mathbf{u}_1|^2 + \frac{1}{2}\alpha_2\rho_2|\mathbf{u}_2|^2 = \frac{1}{2}\rho|\mathbf{u}|^2 + \frac{1}{2}\rho Y_1 Y_2 |\mathbf{W}|^2. \qquad (7) $$

Here, the notations are defined as follows: for $k = 1, 2$, the volume fraction of the $k^{th}$ phase is denoted by $\alpha_k$, while $\rho_k$ denotes its density; the mixture's density is $\rho = \alpha_1\rho_1 + \alpha_2\rho_2$. In this work, we focus on the barotropic setting, the following results are valid in both the isothermal and isentropic case. We consider for each phase an equation of state given by the energy potential $\rho_k \mapsto e_k(\rho_k)$ which depends on its density only. In the isentropic case, the energy potential corresponds to the internal energy while, in the isothermal case, it corresponds to the Helmholtz free energy. The mixture energy potential $e$ is then given by

$$ e(\rho, Y_1, \alpha_1) = Y_1 e_1 \left( \frac{\rho Y_1}{\alpha_1} \right) + Y_2 e_2 \left( \frac{\rho Y_2}{\alpha_2} \right). \qquad (8) $$

The physical Lagrangian is then obtained as the difference between the kinetic and potential energies and thus writes

$$ \mathcal{L} = \frac{1}{2}\rho|\mathbf{u}|^2 + \frac{1}{2}\rho Y_1 Y_2 |\mathbf{W}|^2 - \rho e(\rho, Y_1, \alpha_1). \qquad (9) $$

In addition to the Lagrangian, the constraints, which define the variations associated to each variable, must be specified in order to apply the Stationary Action Principle. Imposing the transport of volume fraction at the barycentric velocity, we set $\alpha_1 \in \mathcal{B}$ according to the previously presented framework. The remaining constraints are the mass conservation for each phase. However, the choice of considering a single trajectory poses some difficulties in imposing these constraints. Indeed, by considering virtual trajectories associated with the barycentric velocity, only the conservation of the total mass can be expressed through these trajectories, but not that of partial masses. In order to overcome this issue, we consider the remaining variables as free variables, i.e. $Y_1, \mathbf{W} \in \Xi$, and we impose the correct continuity equation of phase 1 as a constraint through a Lagrange multiplier $\xi$ which leads to the following contribution to the Action

$$ \int_t \int_{\mathcal{V}(t)} \xi (\partial_t (\rho Y_1) + \nabla \cdot (\rho Y_1 \mathbf{u} + \rho Y_1 Y_2 \mathbf{W})) d\mathbf{x}dt. $$

A simple integration by parts leads to

$$ - \int_t \int_{\mathcal{V}(t)} \rho Y_1 (D_t \xi + Y_2 \mathbf{W} \cdot \nabla \xi) d\mathbf{x}dt + \text{Boundary Terms.} $$

Here, the material derivative $D_t \bullet$ is taken at the barycentric velocity. Since, in the Stationary Action Principle, all variations are assumed to vanish at the boundary [13], any space-time boundary terms in the Action do not contribute to the dynamics of the system. As such, they may be dropped so as to consider the following Lagrangian density

$$ \mathcal{L} = \frac{1}{2}\rho|\mathbf{u}|^2 + \frac{1}{2}\rho Y_1 Y_2 |\mathbf{W}|^2 - \rho e(\rho, Y_1, \alpha_1) - \rho Y_1 (D_t \xi + Y_2 \mathbf{W} \cdot \nabla \xi). \qquad (10) $$
With reference to the generic Stationary Action Principle framework introduced in [13] and previously discussed, we have $\alpha_1 \in \mathcal{B}$ and $Y_1, W, \xi \in \Xi$, and the resulting system writes

$$D_t \alpha_1 = 0, \qquad (11a)$$

$$\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, \qquad (11b)$$

$$\partial_t (\rho Y_1) + \nabla \cdot (\rho Y_1 \mathbf{u}) + \nabla \cdot (\rho Y_1 Y_2 W) = 0, \qquad (11c)$$

$$\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla \bar{p} + \nabla \cdot (\rho Y_1 Y_2 \nabla \xi \otimes W) = 0, \qquad (11d)$$

$$\nabla \xi - W = 0, \qquad (11e)$$

$$D_t \xi + (Y_2 - Y_1) W \cdot \nabla \xi - \frac{1}{2} (Y_2 - Y_1) |W|^2 + h_1 - h_2 = 0, \qquad (11f)$$

where (11a) and (11b) are the constraints, while (11c)-(11f) are the Euler-Lagrange equations. Here, the mixture pressure is $\bar{p} = \alpha_1 p_1 + \alpha_2 p_2$ and $h_k = e_k + p_k/\rho_k$ represents the specific enthalpy — or free enthalpy in the isothermal case — of phase $k = 1, 2$. One can easily verify that (11b) and (11c) are identical to

$$\partial_t (\alpha_k \rho_k) + \nabla \cdot (\alpha_k \rho_k \mathbf{u}_k) = 0, \quad k = 1, 2. \qquad (12)$$

Since $\xi$ is the Lagrange multiplier used to impose the phasic continuity equations, it has no particular physical meaning. As such, we discard it by using (11e) and replacing its evolution equation by an evolution equation for the relative velocity. We do so by taking the material derivative of (11e) and using (11f), which results in the following system

$$D_t \alpha_1 = 0, \qquad (13a)$$

$$\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, \qquad (13b)$$

$$\partial_t (\rho Y_1) + \nabla \cdot (\rho Y_1 \mathbf{u}) + \nabla \cdot (\rho Y_1 Y_2 W) = 0, \qquad (13c)$$

$$\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla \bar{p} + \nabla \cdot (\rho Y_1 Y_2 W \otimes W) = 0, \qquad (13d)$$

$$D_t W + (\nabla \mathbf{u}) W + \nabla (h_1 - h_2 + \frac{1}{2} (Y_2 - Y_1) |W|^2) = 0, \qquad (13e)$$

$$\nabla \times W = 0. \qquad (13f)$$

The relative velocity equation can also be written as

$$\partial_t W + \nabla (H_1 - H_2) + (\partial_x W - (\partial_x W)^T) \mathbf{u} = 0, \qquad (14)$$

with $H_k = h_k + |\mathbf{u}_k|^2/2$. Notice that, since $\nabla \times W = 0$, the last term of (14) is zero, however we keep this term as it ensures the Galilean invariance principle (see also Section 1.3). The supplementary evolution equation of the Hamiltonian corresponds to the classical conservation of total energy, i.e.

$$\partial_t (\alpha_1 \rho_1 E_1 + \alpha_2 \rho_2 E_2) + \nabla \cdot (\alpha_1 \rho_1 H_1 \mathbf{u}_1 + \alpha_2 \rho_2 H_2 \mathbf{u}_2) = 0, \qquad (15)$$

where $E_k = e_k + \frac{1}{2} \mathbf{u}_k^2$ is the total energy per unit of mass of phase $k = 1, 2$. The resulting model (13) has been previously derived through the theory of thermodynamically compatible systems [6, 10, 22, 28] in [24]. Moreover, it was noted that in 1D the model coincides with the Baer-Nunziato model for a specific choice of interfacial closures. Indeed, rewriting system (13) in terms of the phasic velocities $\mathbf{u}_1$ and $\mathbf{u}_2$ instead of $\mathbf{u}$ and $W$, one
obtains

$$ \partial_t \alpha_1 + \mathbf{u}_I \cdot \nabla \alpha_1 = 0 \qquad (16a) $$

$$ \partial_t (\alpha_1 \rho_1) + \nabla \cdot (\alpha_1 \rho_1 \mathbf{u}_1) = 0 \qquad (16b) $$

$$ \partial_t (\alpha_2 \rho_2) + \nabla \cdot (\alpha_2 \rho_2 \mathbf{u}_2) = 0 \qquad (16c) $$

$$ \partial_t (\alpha_1 \rho_1 \mathbf{u}_1) + \nabla \cdot (\alpha_1 \rho_1 \mathbf{u}_1 \otimes \mathbf{u}_1 + \alpha_1 p_1 \mathbf{Id}) - p_I \nabla \alpha_1 = +\rho Y_1 Y_2 [Y_2 \nabla \times \mathbf{u}_1 + Y_1 \nabla \times \mathbf{u}_2] \times \mathbf{W} \qquad (16d) $$

$$ \partial_t (\alpha_2 \rho_2 \mathbf{u}_2) + \nabla \cdot (\alpha_2 \rho_2 \mathbf{u}_2 \otimes \mathbf{u}_2 + \alpha_2 p_2 \mathbf{Id}) + p_I \nabla \alpha_1 = -\rho Y_1 Y_2 [Y_2 \nabla \times \mathbf{u}_1 + Y_1 \nabla \times \mathbf{u}_2] \times \mathbf{W}, \qquad (16e) $$

with $\mathbf{u}_I = \mathbf{u} = Y_1 \mathbf{u}_1 + Y_2 \mathbf{u}_2$ and $p_I = Y_2 p_1 + Y_1 p_2$. As such, in 1D, the derivation based on the Stationary Action Principle yields system (16) which corresponds to a fully closed barotropic Baer-Nunziato model [1]. Moreover, quite remarkably, with the closures of interface quantities $\mathbf{u}_I$ and $p_I$ determined via the variational principle, the non-conservative products are well defined for weak solutions [7]. For multi-dimensional flows, additional “lift” forces are present in the phasic momentum equations. The appearance of lift forces in multi-velocity models obtained by variational principles is also discussed in [8]. The impact of these forces on the wave structure of the model will be discussed in Section 1.3.

As in the variational principle a single virtual trajectory has been used, the model relies on a single reference frame and the relative velocity is an internal variable of the system that contributes, with the thermodynamic variables, to describe the local state of the two-phase mixture. This approach connects with previous modelling approaches based on extended thermodynamics. For instance, in [18, 19], an extended Gibbs potential $g$ is introduced for the mixture

$$ g_{\text{ext}} = Y_1 g_1 (p, T) + Y_2 g_2 (p, T) + \frac{1}{2} Y_1 Y_2 |\mathbf{W}|^2. \qquad (17) $$

For flows which are not necessarily at pressure equilibrium, such as those considered here, the corresponding potential is

$$ e_{\text{ext}} = Y_1 e_1 + Y_2 e_2 + \frac{1}{2} Y_1 Y_2 |\mathbf{W}|^2. \qquad (18) $$

In [24–26], the authors also consider such an extended internal energy per unit of mass $e_{\text{ext}}$. Note however, that in our derivation, the thermodynamic energy and the relative velocity energy have opposite signs in the Lagrangian. This is to ensure that we recover the correct expression of $e_{\text{ext}}$ in the Hamiltonian. Even though the potential energy is often negatively signed in the Lagrangian, when the energies are quadratic, a change of sign may be necessary in order to obtain the physical Hamiltonian, in particular when constraints are imposed through Lagrange multipliers — see for instance [13] for an example of two derivations through the Stationary Action Principle of the same system, but with a different sign for a quadratic energy term in the Lagrangian.

## 1.3. Analysis of eigenvalues and eigenvectors in the multidimensional case

As the one-dimensional version of system (13) has been extensively studied — see for instance [33] — we briefly recall its hyperbolic structure and discuss the effect on the system’s hyperbolicity of the lift forces in the multi-dimensional setting to further the comparison with the Baer-Nunziato models. In the one-dimensional setting, system (13) is hyperbolic with simple eigenvalues

$$ \lambda_1 = u_2 - c_2, \quad \lambda_2 = u_1 - c_1, \quad \lambda_3 = u, \quad \lambda_4 = u_1 + c_1, \quad \lambda_5 = u_2 + c_2, $$

where $c_k^2 = \rho_k \frac{\partial e_k}{\partial \rho_k}$. This directly results from the fact that in the one-dimensional setting, the model can be recast as the barotropic Baer-nunziato model through a change of variables as thus they share the same eigenstructure. As such, the convective eigenvalue $\lambda_3$ is linearly degenerate, while the acoustic eigenvalues $\lambda_{1,2,4,5}$ are genuinely non-linear provided that $\frac{\partial \rho_k c_k}{\partial \rho_k} \neq 0$. As discussed in [33], this is related to the so-called fundamental thermodynamic derivative and, for admissible equations of state, $\frac{\partial \rho_k c_k}{\partial \rho_k} \neq 0$ holds [26].
In the three-dimensional case, the hyperbolicity of both the Baer-Nunziato model and model (13) may be studied assuming slab symmetry and retaining only the x-direction, owing to the models' rotational invariance. For the Baer-Nunziato model, the tangential velocity components of each phase are transported at their own velocity, leading to two additional double eigenvalues, $\lambda = u_1$ and $\lambda = u_2$, and the model remains hyperbolic. However, system (13) differs from the Baer-Nunziato model through the presence of lift forces. These lift forces lead to an exchange between the phasic momenta thus affecting the wave structure of the model and, in particular, the tangential velocity components. Instead of being advected at their own velocity, the tangential components of both phases are now associated to convective eigenvalue $\lambda = u$. For the model to be hyperbolic in the three-dimensional case, its multiplicity should be 5, however, the system lacks two eigenvectors and thus model (13) is only weakly hyperbolic in the three-dimensional setting.

As discussed in [5], weak hyperbolicity appears in several models where an involution constraint, such as (13f), has to be satisfied. In order to restore the strict hyperbolicity of system (13), one may modify the momentum equation as follows

$$ \partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla \bar{p} + \nabla \cdot (\rho Y_1 Y_2 \mathbf{W} \otimes \mathbf{W}) - \rho Y_1 Y_2 (\nabla \times \mathbf{W}) \times \mathbf{W} = \mathbf{0}. \quad (19) $$

Since $\nabla \times \mathbf{W} = \mathbf{0}$ (13f), equation (19) is equivalent to (13d). When considering (19) instead of (13d), system (13)'s eigenvalues are unchanged, but the model is now strictly hyperbolic with the multiple eigenvalue, $\lambda = u$, still being linearly degenerate, and the remaining simple eigenvalues remaining genuinely non-linear under the same thermodynamic condition.

As already mentioned, system (13) is presented in several contributions [24–26] and its derivation in the aforementioned works is justified through the theory of Symmetric Hyperbolic Thermodynamically Compatible systems. An alternative approach is presented in [26], in order to obtain hyperbolicity. The idea is to consider the evolution equation for the relative velocity $\mathbf{W}$ written as

$$ \partial_t \mathbf{W} + \nabla (H_1 - H_2) = -\varepsilon \mathbf{u} \omega, \quad (20) $$

where $\varepsilon$ denotes the pseudo-unity tensor or Levi-Civita symbol and $\omega = \nabla \wedge \mathbf{W}$. As such, the purely multi-dimensional term of the velocity equation is treated as a source term. To do so, the following supplementary equation for $\omega$ has to be satisfied [25, 26]

$$ \partial_t \omega + \nabla \cdot (\omega \otimes \mathbf{u} - \mathbf{u} \otimes \omega) = \mathbf{0}. \quad (21) $$

This choice is justified in [26] with the aim of obtaining a conservative system. The resulting system is strictly hyperbolic; however, in the three-dimensional case, it is characterized by two zero eigenvalues and its hyperbolic part is therefore not frame-invariant.

## 1.4. Discussion concerning the extension to non-barotropic flows

The present work focuses on the barotropic case as the effects related to the presence of several velocities are of mechanical nature rather than of thermodynamic nature. However, the question of its extension to non-barotropic flows naturally arises. As it turns out, such an extension is rather challenging as the present formalism, as well as the other approaches used to derive system (13), run into some difficulties. The aim of the following discussion is to analyse why the modelling frameworks developed so far cannot provide satisfactory results.

As the Stationary Action Principle deals with reversible processes, no entropy is generated inside the fluid particles. Moreover, entropy exchange through heat conduction is irreversible and also creates entropy, thus we assume no entropy exchange between fluid particles<fnref n="2" />. Under these assumptions, the entropy carried by each

<footnote n="2">²Hyperbolic models for heat conduction have been developed. In these models, heat conduction is related to an additional internal variable and has an irreversible character in the absence of dissipative source terms. In this setting, reversible entropy exchange between fluid particles occurs. However, as such a description is not considered here, it is not relevant for the purpose of our work.</footnote>
fluid particle remains constant along its trajectory. In Eulerian coordinates, this translates into

$$ \partial_t (\alpha_k \rho_k s_k) + \nabla \cdot (\alpha_k \rho_k s_k \mathbf{u}_k) = 0. \qquad (22) $$

In the single-velocity case, when following a material control volume that moves with the flow, by definition, no fluid particles cross the control volume's boundary and thus the amount of mass and entropy in the control volume is constant. When relative velocity is taken into account, the situation is quite different as the fluid particles of each phase follow their own velocity. In any single reference frame, such as the one associated to the barycentric velocity $\mathbf{u}$ that we consider, fluid particles of at least one of the phases enter and exit the corresponding control volume. The total amount of entropy, $\rho s = \alpha_1 \rho_1 s_1 + \alpha_2 \rho_2 s_2$ in the control volume then changes over time. In Eulerian coordinates the entropy equations writes

$$ \partial_t (\rho s) + \nabla \cdot (\rho s \mathbf{u} + \rho Y_1 Y_2 (s_1 - s_2) \mathbf{W}) = 0. \qquad (23) $$

The additional entropy flux $\rho Y_1 Y_2 (s_1 - s_2) \mathbf{W}$ is related to the variation of fluid particles entering or exiting the control volume in the barycentric reference frame.

Several extension to non-barotropic flows of model (13) have been proposed in the Symmetric Hyperbolic Thermodynamically Compatible systems framework [24–26]. Among the extensions, the single-temperature [16, 27] and multi-temperature [25, 26] cases have both been investigated. However, it seems that currently, the framework does not allow to impose an equation such as (23) and instead, conservation of the total entropy in the barycentric reference frame is imposed, i.e. the additional entropy flux related to the relative velocity is not taken into account. When doing so, fluid particles which enter the control volume must exchange some entropy with those leaving it in order to satisfy this conservation law. As such, it induces an artificial heat exchange even when both phases are assumed to be in thermal equilibrium with each other, i.e. the phasic temperatures satisfy $T_1 = T_2 =: T$. This artificial heat flux translates in an extra term in the total energy equation which writes

$$ \partial_t (\alpha_1 \rho_1 E_1 + \alpha_2 \rho_2 E_2) + \nabla \cdot (\alpha_1 \rho_1 H_1 \mathbf{u}_1 + \alpha_2 \rho_2 H_2 \mathbf{u}_2 - \rho Y_1 Y_2 (s_1 - s_2) T \mathbf{W}) = 0. \qquad (24) $$

Another specificity of the framework is that non-conservative fluxes proportional to temperature gradients may appear in the phasic momentum equations. As there is no clear physical justification for such fluxes, the extension to non-barotropic flows remains an open question.

Alternatively, in the framework proposed here, one may impose the phasic entropy equations $D_t^k s_k$, $k = 1, 2$ — with $D_t^k$ denoting the material derivative at velocity $\mathbf{u}_k$ — through Lagrange multiplier $\xi_1$ and $\xi_2$. The difficulty is that there is an interaction between the Lagrange multiplier $\xi$ used to impose the phasic mass conservation and the Lagrange multipliers $\xi_1$ and $\xi_2$ and the relative velocity is then a combination of the gradients of all the Lagrange multipliers. Moreover, in the barotropic case, the Stationary Action Principle yields system (11), which we have been able to reduce into system (13), such that the final system only involves physical variables. Unfortunately, to our knowledge, such a reduction is not possible in the non-barotropic setting and calls for an extension of the modelling framework to overcome these difficulties, which is currently under development.

## 2. INCLUSION OF VELOCITY FLUCTUATIONS FOR PSEUDO-TURBULENT FLOWS

When velocity differences are present in the system, pseudo-turbulence is generated and it couples with the phasic velocities. The variational framework based on the Stationary Action Principle provides a natural and elegant way to include such coupling effects ensuring a coherent structure of the model. In order to achieve this goal, we extend the Lagrangian (10) to include the corresponding pseudo-turbulent kinetic energies. First, similarly to classical turbulence modelling, we decompose the phasic velocities in its average $\langle \mathbf{u}_k \rangle$ and a fluctuation $\mathbf{u}_k'$, i.e.

$$ \mathbf{u}_k = \langle \mathbf{u}_k \rangle + \mathbf{u}_k'. \qquad (25) $$
Here, the averaging operator $\langle \bullet \rangle$ is assumed to satisfy the usual properties of Reynolds averaging, and thus $\langle \mathbf{u}'_k \rangle = 0$. As a result, the average phasic specific kinetic energies write

$$ \frac{1}{2} \langle |\mathbf{u}_k|^2 \rangle = \frac{1}{2} |\langle \mathbf{u}_k \rangle|^2 + \frac{1}{2} \langle |(\mathbf{u}'_k)|^2 \rangle = \frac{1}{2} |\langle \mathbf{u}_k \rangle|^2 + k_k, \quad (26) $$

where $k_k$ are the pseudo-turbulent kinetic energies. Next, we consider these supplementary energies to be internal variables that contribute to the thermodynamic description of the two-phase mixture. As such, they are added as potential energies into the previous Lagrangian density (10), which now writes

$$ \mathcal{L} = \frac{1}{2} \rho |\mathbf{u}|^2 + \frac{1}{2} \rho Y_1 Y_2 |\mathbf{W}|^2 - \rho e (\rho, Y_1, \alpha_1) - \rho Y_1 k_1 - \rho Y_2 k_2 - \rho Y_1 (D_t \xi + Y_2 \mathbf{W} \cdot \nabla \xi). \quad (27) $$

Following [30], we assume the transport of the pseudo-turbulent entropies $(\rho Y_k)^{-2/3} k_k$. The transport of turbulent entropies at phasic velocities, and their coupling with the Baer-Nunziato equations, have also been studied in [15]. Here, however, we assume transport at barycentric velocity $\mathbf{u}$ for both turbulent entropies, in accordance with the single reference frame assumption. This leads to imposing the following constraints

$$ D_t ((\rho Y_k)^{-2/3} k_k) = \partial_t ((\rho Y_k)^{-2/3} k_k) + \mathbf{u} \cdot \nabla ((\rho Y_k)^{-2/3} k_k) = 0. \quad (28) $$

These constraints can be explicitly integrated along the virtual trajectories associated to the barycentric velocity, and as such, they are compatible with the previously presented Stationary Action Principle framework, as we simply set $(\rho Y_k)^{-2/3} k_k \in \mathcal{B}$. The application of the generic framework then results in the following system

$$ D_t \alpha_1 = 0, \quad (29a) $$

$$ \partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, \quad (29b) $$

$$ \partial_t (\rho Y_1) + \nabla \cdot (\rho Y_1 \mathbf{u}) + \nabla \cdot (\rho Y_1 Y_2 \mathbf{W}) = 0, \quad (29c) $$

$$ D_t k_1 + \frac{2}{3} k_1 \nabla \cdot \mathbf{u} + \frac{2}{3} \frac{k_1}{\rho Y_1} \nabla \cdot (\rho Y_1 Y_2 \mathbf{W}) = 0, \quad (29d) $$

$$ D_t k_2 + \frac{2}{3} k_2 \nabla \cdot \mathbf{u} - \frac{2}{3} \frac{k_2}{\rho Y_2} \nabla \cdot (\rho Y_1 Y_2 \mathbf{W}) = 0, \quad (29e) $$

$$ \partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla \left( \bar{p} + \frac{2}{3} \rho (Y_1 k_1 + Y_2 k_2) \right) + \nabla \cdot (\rho Y_1 Y_2 \nabla \xi \otimes \mathbf{W}) = 0, \quad (29f) $$

$$ \mathbf{W} - \nabla \xi = 0, \quad (29g) $$

$$ D_t \xi + \frac{1}{2} (Y_2 - Y_1) |\nabla \xi|^2 + (\tilde{h}_1 - \tilde{h}_2) = 0, \quad (29h) $$
where $\tilde{h}_k = h_k + \frac{5}{3}k_k$. Once again, using (29g) and (29h), we discard the Lagrange multiplier $\xi$ and trade its evolution equation for an evolution equation on the relative velocity $\mathbf{W}$. It results in the following system

$$D_t \alpha_1 = 0, \qquad (30a)$$

$$\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, \qquad (30b)$$

$$\partial_t (\rho Y_1) + \nabla \cdot (\rho Y_1 \mathbf{u}) + \nabla \cdot (\rho Y_1 Y_2 \mathbf{W}) = 0, \qquad (30c)$$

$$D_t k_1 + \frac{2}{3} k_1 \nabla \cdot \mathbf{u} + \frac{2}{3} \frac{k_1}{\rho Y_1} \nabla \cdot (\rho Y_1 Y_2 \mathbf{W}) = 0, \qquad (30d)$$

$$D_t k_2 + \frac{2}{3} k_2 \nabla \cdot \mathbf{u} - \frac{2}{3} \frac{k_2}{\rho Y_2} \nabla \cdot (\rho Y_1 Y_2 \mathbf{W}) = 0, \qquad (30e)$$

$$\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla \left( \bar{p} + \frac{2}{3} \rho (Y_1 k_1 + Y_2 k_2) \right) + \nabla \cdot (\rho Y_1 Y_2 \mathbf{W} \otimes \mathbf{W}) = \mathbf{0}, \qquad (30f)$$

$$\partial_t \mathbf{W} + \nabla (\tilde{H}_1 - \tilde{H}_2) + \left( \frac{\partial \mathbf{W}}{\partial \mathbf{x}} - \left( \frac{\partial \mathbf{W}}{\partial \mathbf{x}} \right)^{\top} \right) \mathbf{u} = \mathbf{0}, \qquad (30g)$$

$$\nabla \times \mathbf{W} = \mathbf{0}. \qquad (30h)$$

Here, the terms $\tilde{H}_k = \tilde{h}_k + |\mathbf{u}_k|^2/2$ include the pseudo-turbulent contributions. As such, the system naturally couples the evolution of the relative velocity $\mathbf{W}$ with those of the fluctuation energies $k_k$. System (30) satisfies all the Galilean invariance requirements. In 1D, it is hyperbolic, i.e. it has real eigenvalues and a full set of eigenvectors. The eigenvalues write

$$\lambda_1 = u_2 - \tilde{c}_2, \quad \lambda_2 = u_1 - \tilde{c}_1, \quad \lambda_{3,4,5} = u, \quad \lambda_6 = u_1 + \tilde{c}_1, \quad \lambda_7 = u_2 + \tilde{c}_2. \qquad (31)$$

where $\tilde{c}_k^2 = c_k^2 + \frac{10}{9}k_k$ are the augmented speeds of sound. The acoustic eigenvalues $\lambda_{1,2}$ and $\lambda_{6,7}$ are identical to the acoustic eigenvalues of the turbulent Baer-Nunziato model presented in [15] as the wave structure only differs in its convective part. The characteristic field associated to $\lambda = u$ are linearly degenerate, as such non-conservative products are well defined even through shocks. Finally, the Hamiltonian is given by the total energy

$$\mathcal{H} = \alpha_1 \rho_1 (E_1 + k_1) + \alpha_2 \rho_2 (E_2 + k_2) \qquad (32)$$

and its evolution equation reads as follows

$$\partial_t (\alpha_1 \rho_1 (E_1 + k_1) + \alpha_2 \rho_2 (E_2 + k_2)) + \nabla \cdot (\alpha_1 \rho_1 \tilde{H}_1 \mathbf{u}_1 + \alpha_2 \rho_2 \tilde{H}_2 \mathbf{u}_2) = 0. \qquad (33)$$

In the barotropic context, the total energy — or total free energy depending on whether we consider isentropic or isothermal flows — serves as mathematical entropy. Here, as the mathematical entropy includes the pseudo-turbulent energies, we obtain a system with a coupling between relative velocity and velocity fluctuations for which an inclusive coherent thermodynamic structure is available. In particular, it will allow for the derivation of dissipative source terms in which the energy related to the relative velocity is transferred to the velocity fluctuations. This is the topic of the following section.

### 3. MODELLING DISSIPATIVE EFFECTS

The Stationary Action Principle allows for the derivation of the reversible part of the flow. As such, there are no dissipation effects in the previously derived system. However, the model admits a mathematical entropy given by the Hamiltonian, which allows to add dissipative effects through source terms that are compatible with the second principle of thermodynamics. Here, we focus on the dissipation induced by drag which drives the system towards kinematic equilibrium through viscosity. In the absence of velocity fluctuations, the work
produced by the drag forces can only converted into heat, however, it is reasonable to assume that the drag between the phases induced by the relative velocity generates velocity fluctuations within each phase.

Denoting $R_\phi$ the source term for a generic variable $\phi$, we consider the following equations

$$ \partial_t \mathbf{W} + \nabla (\tilde{H}_1 - \tilde{H}_2) + \left( \frac{\partial \mathbf{W}}{\partial \mathbf{x}} - \left( \frac{\partial \mathbf{W}}{\partial \mathbf{x}} \right)^\top \right) \mathbf{u} = \mathbf{R}_\mathbf{W}, \quad (34a) $$

$$ D_t k_1 + \frac{2}{3} k_1 \nabla \cdot \mathbf{u} + \frac{2}{3} \frac{k_1}{\rho Y_1} \nabla \cdot (\rho Y_1 Y_2 \mathbf{W}) = R_{k_1}, \quad (34b) $$

$$ D_t k_2 + \frac{2}{3} k_2 \nabla \cdot \mathbf{u} - \frac{2}{3} \frac{k_2}{\rho Y_2} \nabla \cdot (\rho Y_1 Y_2 \mathbf{W}) = R_{k_2}, \quad (34c) $$

while keeping the remaining part of system (30) unchanged. The evolution of the Hamiltonian then writes

$$ \partial_t (\alpha_1 \rho_1 (E_1 + k_1) + \alpha_2 \rho_2 (E_2 + k_2)) + \nabla \cdot (\alpha_1 \rho_1 \tilde{H}_1 \mathbf{u}_1 + \alpha_2 \rho_2 \tilde{H}_2 \mathbf{u}_2) = Y_1 Y_2 \mathbf{W} \cdot \mathbf{R}_\mathbf{W} + Y_1 R_{k_1} + Y_2 R_{k_2}. \quad (35) $$

A simple model for drag is a linear relaxation towards zero for the relative velocity with a characteristic relaxation time-scale $\tau > 0$, thus we choose $\mathbf{R}_\mathbf{W}$ as follows

$$ \mathbf{R}_\mathbf{W} = -\frac{1}{\tau} \mathbf{W}. \quad (36) $$

As the corresponding energy is assumed to be entirely transferred to the sub-scale velocity fluctuations, we cancel out the associated mathematical entropy production through the source terms $R_{k_1}$ and $R_{k_2}$. This suggests the following closures

$$ R_{k_1} = \beta_1 \frac{1}{\tau} Y_2 |\mathbf{W}|^2, \quad R_{k_2} = \beta_2 \frac{1}{\tau} Y_1 |\mathbf{W}|^2, \quad (37) $$

where $\beta_1 + \beta_2 = 1$. These coefficients determine which part of the corresponding energy is transferred to each phase and need to be modelled. A simple, yet reasonable, closure is $\beta_1 = \beta_2 = \frac{1}{2}$. Indeed, it leads to the phase which accounts for the minority of the total mass to generate more fluctuations. Since the repartition is based on mass fractions (37), it accounts for both the occupancy rate of the phases, i.e. volume fractions, and their inertia, i.e. density ratios. For instance, in bubbly flows, a greater proportion of the pseudo-turbulent kinetic energy is typically transferred to the bubbles. In contrast, for dispersed droplet flows, the high density ratio generally leads to the generation of pseudo-turbulence in the surrounding gas phase. However, if the droplet volume fraction is sufficiently low, i.e. when only a few droplets are present, the gaseous phase remains on average largely unaffected, resulting in minimal velocity fluctuations.

At equilibrium, the relative velocity between the phases is zero. The previous system with the corresponding relaxation source terms allows for an out-of-equilibrium description. However, for flows near equilibrium, a closed expression of the relative velocity in terms of the other flow variables can be obtained. To do so, we consider an asymptotic expansion in which flows variables are developed in powers of $\tau$, i.e. for a generic flow variable $\phi$, we consider the following formal expansion $\phi = \phi^{(0)} + \tau\phi^{(1)} + \dots$. Injecting this in the relative velocity equation (34a) with relaxation (36) yields

$$ \mathbf{W} = \tau \nabla (\tilde{h}_1 - \tilde{h}_2) + \mathcal{O}(\tau^2). \quad (38) $$
Neglecting all $O(\tau^2)$ terms, the resulting system reads as follows

$$D_t \alpha_1 = 0, \qquad (39a)$$

$$\partial_t \rho + \nabla \cdot (\rho \mathbf{u}) = 0, \qquad (39b)$$

$$\partial_t (\rho Y_1) + \nabla \cdot (\rho Y_1 \mathbf{u}) - \nabla \cdot (\tau \rho Y_1 Y_2 \nabla (\tilde{h}_1 - \tilde{h}_2)) = 0, \qquad (39c)$$

$$\partial_t (\rho \mathbf{u}) + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) + \nabla \cdot \left( \bar{p} + \frac{2}{3} \rho (Y_1 k_1 + Y_2 k_2) \right) = \mathbf{0}, \qquad (39d)$$

$$D_t k_1 + \frac{2}{3} k_1 \nabla \cdot \mathbf{u} - \frac{2}{3} \frac{k_1}{\rho Y_1} \nabla \cdot \left( \tau \rho Y_1 Y_2 \nabla (\tilde{h}_1 - \tilde{h}_2) \right) = \frac{\tau}{2} Y_2 \left| \nabla (\tilde{h}_1 - \tilde{h}_2) \right|^2, \qquad (39e)$$

$$D_t k_2 + \frac{2}{3} k_2 \nabla \cdot \mathbf{u} + \frac{2}{3} \frac{k_1}{\rho Y_2} \nabla \cdot \left( \tau \rho Y_1 Y_2 \nabla (\tilde{h}_1 - \tilde{h}_2) \right) = \frac{\tau}{2} Y_1 \left| \nabla (\tilde{h}_1 - \tilde{h}_2) \right|^2. \qquad (39f)$$

System (39) represents a second-order two-velocity multi-fluid model. As the relative velocity is only present in the second-order fluxes, its first-order part corresponds to the extension to pseudo-turbulence of the barotropic single-velocity, two-pressure multi-fluid models [4]. As such, its first-order part is hyperbolic and has eigenvalues

$$\lambda_1 = u - \tilde{c}_f, \quad \lambda_{2,3,4,5} = u, \quad \lambda_6 = u + \tilde{c}_f, \qquad (40)$$

where $\tilde{c}_f = Y_1 \tilde{c}_1 + Y_2 \tilde{c}_2$ is the extension to pseudo-turbulence of the so-called frozen speed of sound $c_f = Y_1 c_1^2 + Y_2 c_2^2$. One of the difficulties related to multi-velocity models is their loss of hyperbolicity in the single-pressure limit, see [14, 20, 23, 31]. As the velocity non-equilibrium effects are second-order fluxes in (39), the situation is slightly different. Indeed, assuming pressure equilibrium, the expression of the relative velocity at first order in power of series becomes

$$\mathbf{W} = \tau \nabla (\tilde{h}_1 - \tilde{h}_2) = \tau \left( \frac{\alpha_1 - Y_1}{\rho Y_1 Y_2} \nabla p + \frac{5}{3} \nabla (k_1 - k_2) \right). \qquad (41)$$

The limit of instantaneous pressure relaxation obtained through a source term $(p_1 - p_2)/\epsilon$ with $\epsilon \to 0$ in the volume fraction equation (39a), leads to a two-velocity, single-pressure model which is identical to (39) but for which the volume fraction equation (39a) is replaced by the algebraic closure $p_1 - p_2$. The first-order system in $\tau$ remains hyperbolic with eigenvalues

$$\lambda_1 = u - \tilde{c}_w, \quad \lambda_{2,3,4} = u, \quad \lambda_5 = u + \tilde{c}_w, \qquad (42)$$

with $\rho \tilde{c}_w^2 = \rho c_w^2 + \frac{10}{9}(\alpha_1 \rho_1 k_1 + \alpha_2 \rho_2 k_2)$, the mixture speed of sound which is a modification of the so-called Wood sound velocity $c_w$ [34] to account for the velocity fluctuations. We recall that $c_w$ is defined through

$$\frac{1}{\rho c_w^2} = \frac{\alpha_1}{\rho_1 c_1^2} + \frac{\alpha_2}{\rho_2 c_2^2}. \qquad (43)$$

Such models have already been derived in the literature through asymptotic expansions starting from the Baer-Nunziato model, see [12, 17]. The expression for the relative velocity obtained through such techniques is identical to the laminar part of (41) — when $k_1 = k_2 = 0$ — as the starting model does not include any velocity fluctuations. The present approach thus allows for an extension to pseudo-turbulence of such models in a novel and coherent way through the use of variational principles.

## CONCLUSIONS

We have presented the derivation of a two-fluid model through the Stationary Action Principle which includes both relative velocity between the phases and velocity fluctuations. These effects are modelled as internal
variables that contribute to the description of the state of the two-phase mixture. A key novelty of our approach is the simultaneous derivation of both relative velocity and velocity fluctuations in a single framework so as to obtain the coupling effects between these variables. First, we have focused on a simplified model including only relative velocity. While the resulting system has already been derived and studied in previous works [24–26, 33] via the theory of thermodynamically compatible systems [6, 10, 11, 28], we have shown in this contribution an alternative and original derivation based on the Stationary Action Principle. In the one-dimensional case, the derived system corresponds to a fully closed barotropic Baer-Nunziato model with specific closures and, with these closures, non-conservative products are well defined and it is possible therefore to define shock curves.

Next, we have extended the Lagrangian density by including the energies associated with the velocity fluctuations $u'_k$, $k = 1, 2$, so as to obtain a model that accounts for pseudo-turbulent effects caused by the relative velocity. Finally, dissipative source terms have been added compatibly with irreversible thermodynamics to model the energy transfer towards velocity fluctuations. In particular, we have considered the dissipation induced by drag that drives the system towards kinematic equilibrium. An asymptotic analysis has been performed which yields a closed expression of the relative velocity for flows near kinematic equilibrium which in the laminar case recovers previously derived closures [12]. Overall, the variational framework provides a novel, consistent, and thermodynamically grounded way to include pseudo-turbulence effects in multi-fluid models.

Among future developments, we aim at including added mass effects, as in [2]. Moreover, the extension of the model to complete thermodynamics, i.e. non-barotropic flows, is still a major challenge. As discussed in Section 1.4, in order to avoid an artificial heat flux term in the total energy equation, one may consider multiple Lagrange multipliers for the phasic entropy equations. However, as discussed, with this approach, the closure equation for the relative velocity $W$ would result in a non-trivial combination of the gradients of the Lagrange multipliers and would lead to major difficulties to interpret the resulting model. An approach to overcome these issues is under development.

## ACKNOWLEDGEMENTS

The authors would like to acknowledge S. Gavrilyuk and I. Peshkov for the discussions during the Seventh Workshop on Compressible Multiphase Flows as well as the organizers of the workshop.

## REFERENCES

[1] M.R. Baer and J.W. Nunziato. A two-phase mixture theory for the deflagration-to-detonation transition (DDT) in reactive granular materials. *International Journal of Multiphase Flow*, 12:861–889, 1986.

[2] V. Boniou, R.O. Fox, J.W. Posey, and R.W. Houim. A multiphase model for fluid–particle flows with added mass and phase change. *Chemical Engineering Journal*, 502:157967, 2024.

[3] C. Burtea, S. Gavrilyuk, and C. Perrin. Hamilton’s principle of stationary action in multiphase flow modeling. (hal-03146159), feb 2021.

[4] G. Chanteperdrix, P. Villedieu, and J.-P. Vila. A Compressible Model for Separated Two-Phase Flows Computations. In *ASME 2002 Joint U.S.-European Fluids Engineering Division Conference*, pages 809–816. American Society of Mechanical Engineers Digital Collection, 2002.

[5] F. Dhaouadi and S. Gavrilyuk. An Eulerian hyperbolic model for heat transfer derived via hamilton’s principle: analytical and numerical study. *Proceedings of the Royal Society A*, 480(2283):20230440, 2024.

[6] D. Ferrari, I. Peshkov, E. Romenski, and M. Dumbser. A unified HTC multiphase model of continuum mechanics. *Journal of Computational Physics*, 521:113553, 2025.

[7] T. Gallouët, J.-M. Hérard, and N. Seguin. Numerical modeling of two-phase flows using the two-fluid two-pressure approach. *Mathematical Models and Methods in Applied Sciences*, 14(05):663–700, 2004.

[8] S. Gavrilyuk. ‘Uncertainty’ principle in two fluid-mechanics. *ESAIM: Proceedings and Surveys*, 69:47–55, 2020.

[9] J. A. Geurst. Variational principles and two-fluid hydrodynamics of bubbly liquid/gas mixtures. *Physica A: Statistical Mechanics and its Applications*, 135(2):455–486, April 1986.

[10] S.K. Godunov. An interesting class of quasilinear systems. In *Dokl. Akad. Nauk SSSR*, volume 139, pages 521–523, 1961.

[11] S.K. Godunov, T.Y. Mikhailova, and E.I. Romenskii. Systems of thermodynamically coordinated laws of conservation invariant under rotations. *Siberian Mathematical Journal*, 37(4):690–705, 1996.
[12] H. Guillard and F. Duval. A Darcy law for the drift velocity in a two-phase flow model. *Journal of Computational Physics*, 224(1):288-313, 2007.

[13] W. Haegeman, S. Kokh, M Massot, and G. Orlando. A generic framework to derive systems of conservation laws with source terms and its application to heat conduction in fluid flows: An alternative to the method of moments in kinetic theory of gases? *ESAIM: Proceedings and Surveys*, 2025. (to appear, preprint available on HAL hal-05091807v2).

[14] M. Hantke, S. Müller, and L. Grabowsky. News on Baer-Nunziato-type model at pressure equilibrium. *Continuum Mechanics and Thermodynamics*, 33:1-22, 2021.

[15] J.-M. Hérard. Numerical Modelling of Turbulent Two Phase Flows Using the Two Fluid Approach. In *17th AIAA Computational Fluid Dynamics Conference*, Toronto, Canada, 2005. AIAA.

[16] G. La Spina, M. de' Michieli Vitturi, and A. B. Clarke. Transient numerical model of magma ascent dynamics: Application to the explosive eruptions at the Soufrière Hills Volcano. *Journal of Volcanology and Geothermal Research*, 336:118-139, April 2017.

[17] M. Labois, H. Guillard, and M. Grandotto. A five-equation dissipative model for the simulation of two-phase flows. *WIT Transactions on Engineering Sciences*, 59:257-265, 2008.

[18] D. Lhuillier. From molecular mixtures to suspensions of particles. *Journal de Physique II*, 5(1):19-36, 1995.

[19] D. Lhuillier. Internal variables and the non-equilibrium thermodynamics of colloidal suspensions. *Journal of Non-Newtonian Fluid Mechanics*, 96(1-2):19-30, 2001.

[20] D. Lhuillier, C.-H. Chang, and T. Theofanous. On the quest for a hyperbolic effective-field model of disperse flows. *Journal of Fluid Mechanics*, 731, 2013.

[21] S. Mirjalili, S. S. Jain, and M. Dodd. Interface-capturing methods for two-phase flows: An overview and recent developments. *Center for Turbulence Research Annual Research Briefs*, 2017(117-135):13, 2017.

[22] I. Peshkov, M. Pavelka, E. Romenski, and M. Grmela. Continuum mechanics and thermodynamics in the Hamilton and the Godunov-type formulations. *Continuum Mechanics and Thermodynamics*, 30:1343-1378, 2018.

[23] J. D. Ramshaw and J. A. Trapp. Characteristics, Stability, and Short-Wavelength Phenomena in Two-Phase Flow Equation Systems. *Nuclear Science and Engineering*, 66(1):93-102, 1978.

[24] E. Romenski, D. Drikakis, and E. Toro. Conservative Models and Numerical Methods for Compressible Two-Phase Flow. *Journal of Scientific Computing*, 42:68-95, 2010.

[25] E. Romenski, A. Resnyansky, and E. Toro. Conservative hyperbolic formulation for compressible two-phase flow with different phase pressures and temperatures. *Quarterly of applied mathematics*, 65(2):259-279, 2007.

[26] E. Romenski and E.F. Toro. Compressible two-phase flows: two-pressure models and numerical methods. *Comput. Fluid Dyn. J*, 13:403-416, 2004.

[27] Evgeniy Romenski, Alexander Belozerov, and Ilya Peshkov. Conservative formulation for compressible multiphase flows. *Quarterly of Applied Mathematics*, 74, May 2014.

[28] E.I. Romensky. Hyperbolic systems of thermodynamically compatible conservation laws in continuum mechanics. *Mathematical and computer modelling*, 28(10):115-130, 1998.

[29] R. Saurel and R. Abgrall. A multiphase Godunov method for compressible multifluid and multiphase flows. *Journal of Computational Physics*, 150(2):425-467, 1999.

[30] Richard Saurel, Sergey Gavrilyuk, and François Renaud. A multiphase model with internal degrees of freedom: Application to shock-bubble interaction. *Journal of Fluid Mechanics*, 495:283-321, November 2003.

[31] J. H. Stuhmiller. The influence of interfacial pressure forces on the character of two-phase flow model equations. *International Journal of Multiphase Flow*, 3:551-560, 1977.

[32] F. Thein. The Shizuta-Kawashima Condition for the Barotropic SHTC Two Fluid Model. In *XVI International Conference on Hyperbolic Problems: Theory, Numerics, Applications*, pages 361-371. Springer, 2022.

[33] F. Thein, E. Romenski, and M. Dumbser. Exact and numerical solutions of the Riemann problem for a conservative model of compressible two-phase flows. *Journal of Scientific Computing*, 93(3):83, 2022.

[34] A. B. Wood. *A Textbook of Sound*. G. Bell and Sons Ltd., London, 1930.