# An all-topology two-fluid model for two-phase flows derived through Hamilton's Stationary Action Principle

Ward Haegeman, Giuseppe Orlando, Samuel Kokh, Marc Massot

► To cite this version:

Ward Haegeman, Giuseppe Orlando, Samuel Kokh, Marc Massot. An all-topology two-fluid model for two-phase flows derived through Hamilton's Stationary Action Principle. 2025. hal-05249139v2

HAL Id: hal-05249139

https://hal.science/hal-05249139v2

Preprint submitted on 26 Sep 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
# An all-topology two-fluid model for two-phase flows derived through Hamilton's Stationary Action Principle

Ward Haegeman¹,², Giuseppe Orlando², Samuel Kokh³ and Marc Massot²

¹DMPE, ONERA, Université Paris-Saclay, F-91123 Palaiseau, France

²CMAP, CNRS, École polytechnique, Institut Polytechnique de Paris, Route de Saclay, 91120 Palaiseau, France

³Maison de la Simulation, UAR 3441, Bâtiment 565 PC 190 - Digiteo, CEA Saclay 91191, Gif-sur-Yvette, France

**Subject areas:** applied mathematics, mathematical modelling, fluid mechanics

**Keywords:** compressible two-phase flows, multi-fluid models, all-topology models, symmetric hyperbolic systems, Hamilton's principle

**Corresponding author:**

Ward Haegeman

ward.haegeman@polytechnique.edu

September 26, 2025

**Abstract**

We present a novel multi-fluid model for compressible two-phase flows. The model is derived through a newly developed Stationary Action Principle framework. It is fully closed and introduces a new interfacial quantity, the interfacial work. The closures for the interfacial quantities are provided by the variational principle. They are physically sound and well-defined for all type of flow topologies. The model is shown to be hyperbolic, symmetrizable, and admits an entropy conservation law. Its non-conservative products yield uniquely defined jump conditions which are provided. As such, it allows for the proper treatment of weak solutions. In the multi-dimensional setting, the model presents lift forces which are discussed. The model constitutes a sound basis for future numerical simulations.
# 1 Introduction

The class of multi-fluid models is widely accepted among the modelling strategies to accurately describe two-phase flows. Moreover, multi-fluid models are viable in a wide range of flow regimes and can offer practical advantages. For instance, multi-fluid models, obtained through the method of moments, allow for an Eulerian description of disperse flows that may include poly-dispersion and small-scale dynamics [1]. They are a cost-effective alternative to the usual Lagrangian models for which the sampling cost increases drastically with the dimension of the phase space. In the context of separated phase flows, multi-fluid models avoid the need for any explicit interface tracking or reconstruction by relying on a diffuse interface representation, and therefore, they are naturally adapted to handle important topology changes. Moreover, recent contributions have proposed two-scale multi-fluid models to describe both the separated and disperse regimes [2,3]. Yet, there is still no universally accepted basis multi-fluid formulation for two-phase flows and the development of new models remains an active research area. In particular, the development of models that are out of velocity equilibrium remains challenging, yet necessary to accurately describe phenomena such as atomization, slip-flow, added-mass, particle inertia, etc. One of the key challenges is the derivation of a model that is well-posed and allows for the correct treatment of shocks, while being physically valid over a sufficiently wide range of flow regimes.

Among the existing two-velocity models, single-pressure models are generally non-hyperbolic [4–7]. Although a hyperbolic single-pressure model has been recently derived through kinetic theory in the context of disperse flows [8], the model lacks the mathematical foundations to properly treat shocks as its jump conditions are not well-defined. Moreover, it lacks a supplementary conservation law for the entropy. Consequently, as an alternative to single-pressure models, many efforts have been devoted to the development of two-pressure models [9–12]. Among these models, the Baer-Nunziato model [10], initially proposed for deflagration-to-detonation transition in granular materials, has been extended into a class of two-phase flow models which has become quite popular [12–18]. Not relying on any equilibrium assumption, the model’s dependent variables are the phasic flow variables and the volume fraction. The volume fraction advection term, as well as the momentum and energy exchange terms, depend on interfacial quantities, namely the interfacial velocity $u_I$, and the interfacial pressure $p_I$, which must be closed as functions of the dependent variables. The closure problem is common to many two-pressure models and several expressions for the interfacial quantities have been proposed [9, 11, 14, 19–21]. A numerical comparison between some of the suggested closures can be found in [22,23]. Among the different possibilities, the original formulation [10] suggests that the interfacial velocity be taken equal to the velocity of one phase, generally the least compressible one, while the interfacial pressure be taken as the pressure of the other phase.

Although practical, such simple closures are adapted to specific types of flows and may not be sufficiently general to treat complex flows, which may involve important topology changes and even phase inversions. Other closures for the interfacial quantities have also been proposed. For instance for the velocity, the volume-weighted average $u_I = \alpha_1 u_1 + \alpha_2 u_2$, or the mass-weighted average $u_I = Y_1 u_1 + Y_2 u_2$, where $u_k$ is the velocity of the $k^{th}$ phase, while $\alpha_k$ and $Y_k$ are its volume fraction and mass fraction respectively. Although for specific regimes, these may not be the most accurate estimations of the true interfacial velocity, they present the advantage of being well-defined over the whole domain regardless of the flow configuration. Therefore, we refer to models with such interfacial quantities as all-topology models.

Another aspect to consider when closing the interfacial quantities is their impact on the mathematical structure of the resulting model. When the interfacial quantities are functions of only the dependent variables — and not their gradients as in [19,21]— the model is hyperbolic under some non-resonance condition. However, as these interfacial quantities define the inter-phase momentum and energy exchanges, they appear in the governing equations as non-conservative products. As a consequence, jump conditions for shocks may not be well-defined
and different numerical schemes may converge to different solutions, as illustrated in [22]. For the Riemann problem to be well-posed, first, the characteristic field associated to the interfacial velocity should be linearly degenerate in order for the jump conditions to be uniquely defined. Second, the system’s total entropy equation should satisfy a conservation law in order to select admissible weak solutions. The closures which satisfy these conditions for the Baer-Nunziato model have been identified in [14]. In order to obtain a linearly degenerate characteristic field, the interfacial velocity must be chosen as either one of the phasic velocities or as the mass-weighted average velocity. Using the mass-weighted average velocity as interfacial velocity has been proposed in several contributions, in combination with an interfacial pressure that is usually given by the mixture pressure $p_I = \alpha_1 p_1 + \alpha_2 p_2$ [12], or a mixture total pressure $p_I = \alpha_1 p_1 + \alpha_2 p_2 + \rho Y_1 Y_2 (u_1 - u_2)^2$ [13] — where $\rho = \alpha_1 \rho_1 + \alpha_2 \rho_2$ is the total density and $\rho_k$, $k = 1, 2$, are the phasic densities — which is the pressure Riemann invariant and therefore continuous through interfaces. However, these expressions for the interfacial pressure do not satisfy the second requirement for a well-posed Riemann problem, since they do not lead to a conservation law for the system’s entropy. As it turns out, when the interfacial velocity is the mass-weighted average velocity, the only expression for the interfacial pressure that leads to a conservative entropy equation is given by

$$p_I = \frac{Y_2 T_2 p_1 + Y_1 T_1 p_2}{Y_2 T_2 + Y_1 T_1}, \qquad (1.1)$$

where $p_k$ and $T_k$ are respectively the pressure and temperature of the kth phase [14]. However, since the interfacial pressure describes a momentum exchange due to an imbalance of forces at the interface, it is related to mechanical effects rather than thermal ones. As such, it is not expected to depend on the phasic temperatures, and to the best of our knowledge, there is no clear physical interpretation for expression (1.1), which is obtained through purely mathematical considerations by enforcing an entropy conservation law.

An alternative all-topology compressible two-fluid description [24–26] has been derived by means of the theory of thermodynamically compatible systems [27–30]. In its barotropic version, the model has been extensively studied [31,32], while an alternative derivation, and its extension to include velocity fluctuations, through Stationary Action Principle can be found in [33]. In the mono-dimensional setting, the barotropic model is equivalent to a barotropic Baer-Nunziato model with closures $u_I = Y_1 u_1 + Y_2 u_2$ and $p_I = Y_2 p_1 + Y_1 p_2$, but expressed in terms of total momentum, $\rho u = \alpha_1 \rho_1 u_1 + \alpha_2 \rho_2 u_2$ — notice that here $u = u_I$ — and the relative velocity, $W = u_1 - u_2$, instead of the phasic momenta $\alpha_k \rho_k u_k$, $k = 1, 2$. In the multi-dimensional setting, compared to the Baer-Nunziato model, additional lift-like forces appear, as well as the involution constraint $\nabla \wedge W = 0$. However, as discussed in [33], the main challenge related to this model is its extension to full thermodynamics, i.e. to non-barotropic flows. Several versions of the model which include full thermodynamics have been proposed, as either single- or two-temperature models [24–26], but all share the same type of underlying assumption, which is the conservation of entropy with mixture velocity $u$. For smooth solutions of the Euler equations, the entropy density is conserved by the flow. Extending this to two-fluid models, one may assume the conservation of the phasic entropy density, $\alpha_k \rho_k s_k$, $k = 1, 2$, following its own velocity $u_k$. This then leads to the following equation for the total entropy density, $\rho s = \alpha_1 \rho_1 s_1 + \alpha_2 \rho_2 s_2$, written in terms of mixture velocity $u$ and relative velocity $W$,

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

In the framework of thermodynamically compatible systems, phasic entropies are assumed to be conserved in the frame associated to the mixture velocity instead of their own velocity. This leads to a total entropy equation for which the entropy flux proportional to $W$ in (1.2) is absent [24–26,34]. Consequently, non-zero relative velocity induces an entropy exchange between the phases which induces an artificial heat flux in the total energy equation, which then writes

$$\partial_t (\alpha_1 \rho_1 E_1 + \alpha_2 \rho_2 E_2) + \nabla \cdot (\alpha_1 \rho_1 H_1 u_1 + \alpha_2 \rho_2 H_2 u_2 - \rho Y_1 Y_2 (s_1 T_1 - s_2 T_2) W) = 0, \qquad (1.3)$$
where $E_k = e_k + \mathbf{u}_k^2/2$, and $H_k = h_k + \mathbf{u}_k^2/2$, are respectively the specific total energy and specific total enthalpy of the $k^{th}$ phase, with $e_k$, and $h_k = e_k + p_k/\rho_k$, its specific internal energy and specific enthalpy. Note that no entropy production is associated to this heat flux term and that it is not related to any heat diffusion phenomenon. Moreover, the framework then also leads to the appearance of non-conservative products proportional to the phasic temperature gradients in the phasic momentum equations and for which there is no clear physical interpretation. As a result, the quest for an *all-topology* two-fluid model, which has both the mathematical properties required to ensure its well-posedness and to properly treat shocks, as well as the physical consistency that allows for a clear interpretation of its various terms, does not seem to have received a conclusive answer yet.

In this contribution, we propose a novel *all-topology* two-fluid model derived through Hamilton's Stationary Action Principle. Compared to the variational principle used in [33], which relies on a single family of trajectories, here the variational principle relies on two family of trajectories, one for each phase. The resulting model is fully closed and the resulting expressions for the interfacial quantities recover previously used and physically sound expressions. Moreover, a new interfacial quantity, representing the inter-phase work produced at the interface, appears, and is coined interfacial work accordingly. The model is analysed in the one-dimensional setting, under a non-resonance condition — which is classical for such two-fluid models [9,14] — the model is hyperbolic and symmetrizable, which ensures its local well-posedness. Furthermore, owing to the linearly degenerate structure of the interfacial wave, its non-conservative products are well-defined for weak solutions and thus the model admits a unique set of jump conditions. The conservation law satisfied by the total entropy allows to select admissible shocks in agreement with the second principle of thermodynamics.

In the multi-dimensional setting, an additional force, akin to lift, appears in the phasic momentum equations and is a direct consequence of the imposed interfacial velocity in the Hamiltonian Principle. Its expression implies the relative velocity between the two phases but also a new potential velocity field which has its own evolution equation, and is generated by the pressure non-equilibrium between the two phases. The presence of lift forces has also been obtained in the thermodynamically compatible two-phase models [24–26], and in the Hamiltonian models found in [35], but are not yet well-understood. We compare their expression in both formulations but the proper modelling of such effects requires further studies.

The outline of this contribution is as follows. We end this section by specifying some notations and vector analysis conventions. Next, in Section 2, the two-trajectory Stationary Action Principle framework is presented. Then, in Section 3, we derive and discuss the new model. Its mathematical analysis in the one-dimensional setting is conducted in Section 4, some comments on the multi-dimensional setting are also provided. Final remarks and conclusions are presented in Section 5. Appendices A and B provide supplementary details concerning the modelling framework.

We adopt the following conventions. Vectors of $\mathbb{R}^d$ are represented in column form, and the identity matrix is denoted $I_d$. The derivative of a scalar field with respect to a vector field is its differential and is represented by a row-vector, its gradient is the transpose of its differential and is a column-vector. Similarly, the derivative of a vector with respect to another vector is its Jacobian matrix. The gradient of a vector field is defined as the transpose of its Jacobian matrix. The divergence of a matrix is the column vector obtained by applying the divergence operator to each row-vector of the matrix.

For a phasic velocity field $\mathbf{u}_k$, the associated material derivative is denoted $D_t^k$, i.e. $D_t^k \bullet = \partial_t \bullet + (\mathbf{u}_k \cdot \nabla)\bullet$. The mass-weighted averaged velocity $Y_1\mathbf{u}_1 + Y_2\mathbf{u}_2$ is denoted $\mathbf{u}$ as indicated previously, it will also be referred to as mixture velocity. The material derivative operator defined by this velocity is simply denoted $D_t$. Interfacial velocities are denoted $\mathbf{u}_I$, when the interfacial velocity equals the mixture velocity, its index may be dropped without ambiguity.
## 2 The Stationary Action Principle with two families of trajectories

Hamilton's Stationary Action Principle allows for a variational derivation of equation of motions under the reversibility assumption. Its application to fluid mechanics allows to derive the equations for ideal fluids in both Lagrangian and Eulerian coordinates [36–38]. The Lagrangian derivations are based on the introduction of fluid particle trajectories $x = \Phi(t, X)$, where $x$ and $X$ are the Eulerian and Lagrangian coordinates respectively. The perturbations of the different flow variables are then related to perturbations of the trajectories and naturally account for constraints such as mass conservation. On the other hand, derivations in the Eulerian setting assume independent variations for all flow variables and constraints such as mass conservation are included by means of Lagrange multipliers. However, additional constraints, called Lin's constraints, must be added in order to obtain the correct equations of motion. As noted in [36], this is related to the fact that the Hamiltonian principle requires to maintain fixed end points during the variational process. In the Lagrangian derivation, this is naturally included by considering compactly supported perturbations of the trajectories. However, in the Eulerian derivation, compactly supported velocity perturbations do not necessarily maintain fixed end points and thus additional constraints are added. It then results in an expression of the velocity field in terms of Clebsch potentials from which the equations of motion must be deduced. In the Lagrangian derivation, one immediately obtains the correct equations of motion, and thus the Lagrangian approach is preferred in this work.

A generic Stationary Action Principle framework has been presented in [39] and used for the derivation of the hyperbolic model for heat transfer presented in [40]. The framework assumes a single family of trajectories and is thus well-adapted for single-velocity models. Extension to models without velocity equilibrium has been achieved in [33] by assuming the single trajectory to be a virtual trajectory associated the mass-weighted average velocity and by adding the relative velocity as an internal variable associated to phasic mass conservation constraints. This has allowed for an original derivation of the barotropic version of the two-phase flow model originally obtained through the theory of thermodynamically compatible systems in [24–26], as well as its extension to include pseudo-turbulence. However, as mentioned in the introduction, the lifting of the barotropic assumption in this single-trajectory framework leads to artificial heat exchanges.

The introduction of two family of trajectories, $\Phi_1$ and $\Phi_2$, has allowed for the application of the variational principle to fluid mixtures in the context of two-phase flows. Single-pressure models, classically obtained through phase averaging, have been obtained by application of the Stationary Action Principle to averaged energy densities [41]. Interactions between the velocity fields such as added-mass effects have also been considered [42]. Two-pressure models have also been derived [43] using one of the phasic velocities to govern the volume fraction dynamics.

In this section, we extend the generic framework presented in [39] into a two-trajectories framework. To do so, we consider a generic Lagrangian density $\mathcal{L}$, which is assumed of the following form,

$$ \mathcal{L} = \mathcal{L}_1(m_1, u_1, b_{1,i}, b_{2,i}, \xi_j, \nabla \xi_j, D_t^1 \xi_j) + \mathcal{L}_2(m_2, u_2, b_{2,i}, b_{1,i}, \xi_j, \nabla \xi_j, D_t^2 \xi_j). \quad (2.1) $$

Here, for $k = 1, 2$, $u_k$ denotes the velocity associated to the trajectory $\Phi_k$, while $m_k$ is the partial mass that is conserved along this trajectory. Variables $b_{k,i}$, $i = 1, \dots, n_k$ are transported at velocity $u_k$, they allow to represent constraints, which may be expressed as transport equations. As such, $b_{k,i}$ may play the role of a phasic entropy, mass fractions in the context of multi-component phases, turbulent entropies [19,33,44], etc. Finally, variables $\xi_j$, $j = 1, \dots, m$, are free variables in the sense that they are not related to any of the trajectories as they will have their own independent variations. These may represent coupling variables, sub-scale phenomena such as micro-intertia [43], or be related to internal constraints. We recall that the notation $D_t^k \bullet$ refers to the material derivative associated to the velocity $u_k$.
The Stationary Action Principle then dictates that the equations of motion are critical points of the Action A, which is the space-time integral of the Lagrangian density L. As such, the generic equations are determined by

$$ \delta A = \delta \int_t \int_{\Omega_t} \mathcal{L} d\mathbf{x} dt = 0, \qquad (2.2) $$

subjected to the following constraints, related to the conservation of the partial masses,

$$ \partial_t m_1 + \nabla \cdot (m_1 \mathbf{u}_1) = 0, \qquad (2.3a) $$

$$ \partial_t m_2 + \nabla \cdot (m_2 \mathbf{u}_2) = 0, \qquad (2.3b) $$

and the transport of the advected quantities,

$$ D_t^1 b_{1,i} = 0, \quad i = 1, \dots, n_1, \qquad (2.3c) $$

$$ D_t^2 b_{2,i} = 0, \quad i = 1, \dots, n_2. \qquad (2.3d) $$

The space of perturbations is then parametrized by a one-parameter family of perturbed trajectories $\Phi_k(t, \mathbf{X}, \varepsilon)$. These then define a virtual displacement, $\delta\Phi_k(t, \mathbf{x})$, that we express in Eulerian coordinates. Perturbations of the velocities $\delta\mathbf{u}_k$, partial masses $\delta m_k$, and advected quantities $\delta b_{k,i}$ are then related to the virtual motion $\delta\Phi_k$, through the following relations,

$$ \delta \mathbf{u}_k = D_t^k (\delta \mathbf{\Phi}_k) - (\nabla \mathbf{u}_k)^\top \delta \mathbf{\Phi}_k, \qquad (2.4a) $$

$$ \delta m_k = -\nabla \cdot (m_k \delta \mathbf{\Phi}_k), \qquad (2.4b) $$

$$ \delta b_{k,i} = -(\delta \mathbf{\Phi}_k)^\top \nabla b_{k,i}. \qquad (2.4c) $$

These perturbations encode the constraints (2.3), and their expressions are classical, see for instance [36,39,40,43].

Variables $\xi_j$ are assumed to have independent variations $\delta\xi_j$. Commuting the $\delta$ operator with space and time derivatives, then yields the following expressions for the remaining variations,

$$ \delta (\nabla \xi_j) = \nabla (\delta \xi_j), \qquad (2.4d) $$

$$ \delta (D_t^k \xi_j) = D_t^k (\delta \xi_j) + (\nabla \xi_j)^\top \delta \mathbf{u}_k. \qquad (2.4e) $$

Finally, following the Hamiltonian principle according to which all boundary points are kept fixed, the perturbations are assumed to be with compact support in the space-time domain. We then obtain the following results, proofs of which are reported in appendix A.

**Theorem 2.1** (The generic system). *The Stationary Action Principle (2.2), applied to the generic Lagrangian (2.1), with variations (2.4), corresponding to the constraints (2.3), yields the following evolution equations,*

$$ \partial_t \mathbf{K}_1 + \nabla \cdot (\mathbf{K}_1 \otimes \mathbf{u}_1 + \mathbf{\Pi}_1) = - \sum_{j=1}^{m} \left[ \partial_t \mathcal{M}_{1,j} + \nabla \cdot (\mathcal{M}_{1,j} \mathbf{u}_1) + \nabla \cdot \mathbf{D}_{1,j} - \frac{\partial \mathcal{L}_1}{\partial \xi_j} \right] \nabla \xi_j \\ + \sum_{i=1}^{n_2} \frac{\partial \mathcal{L}_1}{\partial b_{2,i}} \nabla b_{2,i} - \sum_{i=1}^{n_1} \frac{\partial \mathcal{L}_2}{\partial b_{1,i}} \nabla b_{1,i}, \qquad (2.5a) $$

$$ \partial_t \mathbf{K}_2 + \nabla \cdot (\mathbf{K}_2 \otimes \mathbf{u}_2 + \mathbf{\Pi}_2) = - \sum_{j=1}^{m} \left[ \partial_t \mathcal{M}_{2,j} + \nabla \cdot (\mathcal{M}_{2,j} \mathbf{u}_2) + \nabla \cdot \mathbf{D}_{2,j} - \frac{\partial \mathcal{L}_2}{\partial \xi_j} \right] \nabla \xi_j \\ - \sum_{i=1}^{n_2} \frac{\partial \mathcal{L}_1}{\partial b_{2,i}} \nabla b_{2,i} + \sum_{i=1}^{n_1} \frac{\partial \mathcal{L}_2}{\partial b_{1,i}} \nabla b_{1,i}, \qquad (2.5b) $$
and

$$ \partial_t (\mathcal{M}_{1,j} + \mathcal{M}_{2,j}) + \nabla \cdot (\mathcal{M}_{1,j} \mathbf{u}_1 + \mathcal{M}_{2,j} \mathbf{u}_2 + \mathbf{D}_{1,j} + \mathbf{D}_{2,j}) - \frac{\partial \mathcal{L}_1}{\partial \xi_j} - \frac{\partial \mathcal{L}_2}{\partial \xi_j} = 0, \quad (2.5c) $$

for $j = 1, \dots, m$. Here, $\mathbf{K}_1$ and $\mathbf{K}_2$ are the generalized partial momenta, while $\mathbf{\Pi}_1$ and $\mathbf{\Pi}_2$ are their associated generalized pressure tensors, they are defined as

$$ \mathbf{K}_k = \left( \frac{\partial \mathcal{L}_k}{\partial \mathbf{u}_k} \right)^{\top}, \quad \mathbf{\Pi}_k = \left( \mathcal{L}_k - m_k \frac{\partial \mathcal{L}_k}{\partial m_k} \right) \mathbf{I}_d - \sum_{j=1}^{m} \nabla \xi_j \otimes \mathbf{D}_{k,j}, \quad k = 1, 2, \quad (2.6a) $$

The Euler-Lagrange equations related to the free variables are expressed through their conjugate variables,

$$ \mathcal{M}_{k,j} = \frac{\partial \mathcal{L}_k}{\partial (\mathbf{D}_t^k \xi_j)}, \quad \mathbf{D}_{k,j} = \left( \frac{\partial \mathcal{L}_k}{\partial (\nabla \xi_j)} \right)^{\top}, \quad j = 1, \dots, m, \quad k = 1, 2. \quad (2.6b) $$

**Remark 2.1.** The terms proportional to $\nabla \xi_j$ on the right-hand sides of eq. (2.5a) and eq. (2.5b) sum up to zero owing to eq. (2.5c), the remaining part of the right-hand sides also sum up to zero owing to their anti-symmetry, therefore the generic system ensures the conservation of the total momentum, $\mathbf{K}_1 + \mathbf{K}_2$, in agreement with Noether's theorem.

**Remark 2.2.** When $\mathcal{L}_1$ and $\mathcal{L}_2$ do not depend on $\mathbf{D}_t^k \xi_j$ and $\mathbf{D}_t^2 \xi_j$, eq. (2.5c) yields an algebraic constraint instead of an evolution equation.

For convenience, in the following, we will denote by $\mathcal{B}_1 = \{b_{1,i}, i = 1, \dots, n_1\}$ and $\mathcal{B}_2 = \{b_{2,i}, i = 1, \dots, n_2\}$, the sets of transported variables, while $\Xi = \{\xi_j, j = 1, \dots, m\}$ will denote the set of free variables.

As the generic Lagrangian (2.1) does not explicitly depend the time variable, the system's total energy, also called Hamiltonian, is conserved as specified by the following result.

**Theorem 2.2** (Conservation of the Hamiltonian). For $k = 1, 2$, define the partial Hamiltonian $\mathcal{H}_k$ as follows,

$$ \mathcal{H}_k = \mathbf{K}_k^\top \mathbf{u}_k + \sum_{j=1}^{m} \mathcal{M}_{k,j} \mathbf{D}_t^k \xi_j - \mathcal{L}_k. \quad (2.7) $$

The generic system composed of eqs. (2.3) and eqs. (2.5) admits a supplementary conservation law satisfied by the Hamiltonian, $\mathcal{H} = \mathcal{H}_1 + \mathcal{H}_2$, which writes

$$ \partial_t \mathcal{H} + \nabla \cdot \left( \mathcal{H}_1 \mathbf{u}_1 + \mathcal{H}_2 \mathbf{u}_2 + \mathbf{\Pi}_1^\top \mathbf{u}_1 + \mathbf{\Pi}_2^\top \mathbf{u}_2 + \sum_{j=1}^{m} (\mathbf{D}_{1,j} \mathbf{D}_t^1 \xi_j + \mathbf{D}_{2,j} \mathbf{D}_t^2 \xi_j) \right) = 0. \quad (2.8) $$

The equation for the partial Hamiltonian $\mathcal{H}_k$ writes

$$ \partial_t \mathcal{H}_k + \nabla \cdot \left( \mathcal{H}_k \mathbf{u}_k + \mathbf{\Pi}_k^\top \mathbf{u}_k + \sum_{j=1}^{m} \mathbf{D}_{k,j} \mathbf{D}_t^k \xi_j \right) = \sum_{i=1}^{n_k} \frac{\partial \mathcal{L}_{k,i}}{\partial b_{k,i}} \partial_t b_{k,i} - \sum_{i=1}^{n_{k'}} \frac{\partial \mathcal{L}_k}{\partial b_{k',i}} \partial_t b_{k',i} \\ + \sum_{j=1}^{m} \left[ \partial_t \mathcal{M}_{k,j} + \nabla \cdot (\mathcal{M}_{k,j} \mathbf{u}_k) + \nabla \cdot \mathbf{D}_{k,j} - \frac{\partial \mathcal{L}_k}{\partial \xi_j} \right] \partial_t \xi_j, \quad (2.9) $$

with $k' = 3 - k$, and $k = 1, 2$.
We conclude this section by showing how the present framework may be applied to derive the Baer-Nunziato two-phase flow model with its original closures [10]. We consider two compressible and immiscible fluids, each equipped with a complete equation of state, given by a strictly concave function $(\tau_k, e_k) \mapsto s_k$, with $\tau_k = 1/\rho_k$, the specific volume. The Gibbs identity, $T_k ds_k = de_k + p_k d\tau_k$, then defines the phasic temperature, $T_k > 0$, and the phasic pressure, $p_k$. The $k^{th}$ phase's occupancy rate is given by its volume fraction $\alpha_k \in (0, 1)$, which satisfies the saturation constraint $\alpha_1 + \alpha_2 = 1$. The Lagrangian is given by the difference between the kinetic energy and the potential energy which, in the absence of external forces reduces to the internal thermodynamic energy. Hence, we obtain

$$ \mathcal{L} = \frac{1}{2} m_1 |\mathbf{u}_1|^2 - m_1 e_1 \left( \frac{m_1}{\alpha_1}, s_1 \right) + \frac{1}{2} m_2 |\mathbf{u}_2|^2 - m_2 e_2 \left( \frac{m_2}{1 - \alpha_1}, s_2 \right), \quad (2.10) $$

where $m_k = \alpha_k \rho_k$ is the partial mass of the $k^{th}$ phase. In the absence of heat exchange and under the reversibility assumption, fluid particles do not exchange entropy, nor is there any entropy creation, as such each fluid particle is assumed to carry a fixed amount of entropy. Moreover, following the original closures, we impose the interfacial velocity $\mathbf{u}_I = \mathbf{u}_1$ for the volume fraction transport. As a result, we consider the Lagrangian (2.10) with $\mathcal{B}_1 = \{\alpha_1, s_1\}$, $\mathcal{B}_2 = \{s_2\}$, and $\Xi = \emptyset$. Application of Theorem 2.1 then yields the phasic momentum equations, while application of Theorem 2.2 allows to exchange the phasic entropy equations for the phasic total energy equations as $\mathcal{H}_k = \alpha_k \rho_k E_k$. The complete system that results then writes

$$ \partial_t \alpha_1 + \mathbf{u}_1 \cdot \nabla \alpha_1 = 0, \quad (2.11a) $$

$$ \partial_t (\alpha_1 \rho_1) + \nabla \cdot (\alpha_1 \rho_1 \mathbf{u}_1) = 0, \quad (2.11b) $$

$$ \partial_t (\alpha_2 \rho_2) + \nabla \cdot (\alpha_2 \rho_2 \mathbf{u}_2) = 0, \quad (2.11c) $$

$$ \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{I}_d) - p_2 \nabla \alpha_1 = 0, \quad (2.11d) $$

$$ \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{I}_d) + p_2 \nabla \alpha_1 = 0, \quad (2.11e) $$

$$ \partial_t (\alpha_1 \rho_1 E_1) + \nabla \cdot (\alpha_1 \rho_1 H_1 \mathbf{u}_1) - p_2 \mathbf{u}_1 \cdot \nabla \alpha_1 = 0, \quad (2.11f) $$

$$ \partial_t (\alpha_2 \rho_2 E_2) + \nabla \cdot (\alpha_2 \rho_2 H_2 \mathbf{u}_2) + p_2 \mathbf{u}_1 \cdot \nabla \alpha_1 = 0. \quad (2.11g) $$

This is exactly the Baer-Nunziato system [10] with interfacial velocity $\mathbf{u}_I = \mathbf{u}_1$ and interfacial pressure $p_I = p_2$. It is important to note that here, only the choice of $\mathbf{u}_1$ as interfacial velocity is imposed as a constraint, since $\alpha_1 \in \mathcal{B}_1$, however, the corresponding closure for the interfacial pressure results directly from the variational principle.

Much of the physical information about the system is contained within the Lagrangian, yet, the constraints associated to the variables play an equally important role. For instance, if we change the volume fraction constraint and set it to be a free variable, i.e. $\alpha_1 \in \Xi$, instead of eq. (2.11a), we obtain an algebraic constraint in accordance with Remark 2.2. This constraint reads $p_1 = p_2$, and results in a single-pressure model that is non-hyperbolic.

We now move on to the next section in which an all-topology two-fluid model is derived.

# 3 An all-topology two-phase flow model

## 3.1 Derivation of the model

We consider the same physical setting and assumptions as used in the previous section to derive the Baer-Nunziato model through Stationary Action Principle. In order to derive an all-topology two-fluid model, the constraint on the volume fraction must be modified as it leads to the interfacial velocity to be equal to one of the phasic velocities.

Following [12,13], we consider the mass-weighted average velocity, $\mathbf{u} = Y_1 \mathbf{u}_1 + Y_2 \mathbf{u}_2$, as interfacial velocity. When drag is considered, this is the velocity to which the phasic velocities relax [12]. It is also the choice that is made in the two-phase flow models derived by the theory
of thermodynamically compatible systems [24-26]. Moreover, as shown in [14], under mild assumptions, this is the only all-topology velocity closure which allows the non-conservative products in the Baer-Nunziato model to be uniquely defined in the case of discontinuities. Therefore, there are both good physical and mathematical arguments in favour of this closure.

However, this new constraint related to the volume fraction cannot be exactly integrated along the phasic trajectories, and as such, no exact expression for the perturbation $\delta\alpha_1$ as a function of $\delta\Phi_1$ and $\delta\Phi_2$ is available. To circumvent this issue, we impose this new constraint through a Lagrange multiplier $\zeta$. This leads to a supplementary contribution, $\zeta\rho D_t\alpha_1$, to the Lagrangian, which we write as

$$ \zeta\rho D_t\alpha_1 = \zeta m_1 D_t^1\alpha_1 + \zeta m_2 D_t^2\alpha_1. \quad (3.1) $$

Adding this contribution to the physical Lagrangian (2.10), results in the following effective Lagrangian,

$$ \mathcal{L} = \frac{1}{2}m_1 |\mathbf{u}_1|^2 - m_1 e_1 \left( \frac{m_1}{\alpha_1}, s_1 \right) + \zeta m_1 D_t^1 \alpha_1 \\ + \frac{1}{2}m_2 |\mathbf{u}_2|^2 - m_2 e_2 \left( \frac{m_2}{1-\alpha_1}, s_2 \right) + \zeta m_2 D_t^2 \alpha_1. \quad (3.2) $$

The previously proposed framework is then applied with $\mathcal{L}_1$ and $\mathcal{L}_2$ being given by the first and second lines of eq. (3.2)'s right hand-side respectively and where the constraints are given by $\Xi = \{\alpha_1, \zeta\}$ and $\mathcal{B}_k = \{s_k\}, k = 1, 2$.

Application of Theorem 2.1 then yields the following equations,

$$ \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{I}_d) - p_I \nabla \alpha_1 = +\mathbf{f}, \quad (3.3a) $$

$$ \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{I}_d) + p_I \nabla \alpha_1 = -\mathbf{f}, \quad (3.3b) $$

$$ \partial_t \alpha_1 + \mathbf{u} \cdot \nabla \alpha_1 = 0, \quad (3.3c) $$

$$ \partial_t \zeta + \mathbf{u} \cdot \nabla \zeta - \frac{p_1 - p_2}{\rho} = 0, \quad (3.3d) $$

with the following expressions for the interfacial pressure $p_I$, and the additional momentum exchange force $\mathbf{f}$,

$$ p_I = Y_2 p_1 + Y_1 p_2, \quad (3.4a) $$

$$ \mathbf{f} = \rho Y_1 Y_2 [\nabla \zeta \otimes \nabla \alpha_1 - \nabla \alpha_1 \otimes \nabla \zeta] \mathbf{W}, \quad (3.4b) $$

where we recall that $\mathbf{W} = \mathbf{u}_1 - \mathbf{u}_2$ denotes the relative velocity. These equations are to be supplemented with

$$ \partial_t (\alpha_1 \rho_1) + \nabla \cdot (\alpha_1 \rho_1 \mathbf{u}_1) = 0, \quad (3.5a) $$

$$ \partial_t (\alpha_2 \rho_2) + \nabla \cdot (\alpha_2 \rho_2 \mathbf{u}_2) = 0, \quad (3.5b) $$

$$ \partial_t s_1 + \mathbf{u}_1 \cdot \nabla s_1 = 0, \quad (3.5c) $$

$$ \partial_t s_2 + \mathbf{u}_2 \cdot \nabla s_2 = 0, \quad (3.5d) $$

to form a complete and closed system. This systems admits a supplementary conservation law for the Hamiltonian, however, as is custom, we replace the phasic entropy equations by the phasic energy equations, as given by Theorem 2.2, and use the total entropy balance as supplementary equation. Moreover, notice that the expression of $\mathbf{f}$ involves two gradients, $\nabla\alpha_1$ and $\nabla\zeta$, and so the system, as is, is not in the form of a first-order system of balance laws. To recast it in such a form, we introduce $\mathbf{v} = \nabla\zeta$ as a new independent variable which replaces $\zeta$. The resulting
first-order system then writes

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

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

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

$$ \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{I}_d) - p_I \nabla \alpha_1 = +\mathbf{f}, \qquad (3.6d) $$

$$ \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{I}_d) + p_I \nabla \alpha_1 = -\mathbf{f}, \qquad (3.6e) $$

$$ \partial_t (\alpha_1 \rho_1 E_1) + \nabla \cdot (\alpha_1 \rho_1 H_1) - (p\mathbf{u})_I \cdot \nabla \alpha_1 = +\mathbf{u} \cdot \mathbf{f}, \qquad (3.6f) $$

$$ \partial_t (\alpha_2 \rho_2 E_2) + \nabla \cdot (\alpha_2 \rho_2 H_2) + (p\mathbf{u})_I \cdot \nabla \alpha_1 = -\mathbf{u} \cdot \mathbf{f}, \qquad (3.6g) $$

$$ \partial_t \mathbf{v} + (\mathbf{u} \cdot \nabla) \mathbf{v} + (\nabla \mathbf{u}) \mathbf{v} - \nabla \left( \frac{p_1 - p_2}{\rho} \right) = 0, \qquad (3.6h) $$

with the new interfacial quantity, the interfacial work $(p\mathbf{u})_I$, defined by

$$ (p\mathbf{u})_I = Y_1 p_2 \mathbf{u}_1 + Y_2 p_1 \mathbf{u}_2, \qquad (3.7) $$

and the right-hand side that is defined by

$$ \mathbf{f} = \rho Y_1 Y_2 [\mathbf{v} \otimes \nabla \alpha_1 - \nabla \alpha_1 \otimes \mathbf{v}] \mathbf{W}. \qquad (3.8) $$

This system must then be adjoined with the condition $\mathbf{rot}(\mathbf{v}) = \mathbf{0}$. As this constraint is preserved in time by (3.6h), it reduces to a condition on the initial condition for $\mathbf{v}$. Moreover, one can show that the system is Galilean invariant.

Before moving on to the discussion about the obtained model's features, we mention that a fully Eulerian derivation is also possible. To do so, one considers the following Lagrangian,

$$ \mathcal{L} = \sum_{k=1}^{2} \left\{ \frac{1}{2} m_k \mathbf{u}_k^2 - m_k e_k \left( \frac{m_k}{\alpha_k}, s_k \right) + \zeta m_k D_t^k \alpha_1 \right. \\ \left. + \xi_k \left[ \partial_t m_k + \nabla \cdot (m_k \mathbf{u}_k) \right] + m_k \theta_k D_t^k s_k + m_k \sum_i \lambda_{k,i} D_t^k \mu_{k,i} \right\}. \qquad (3.9) $$

Its first line corresponds to the effective Lagrangian (3.2), whereas the second line contains the phasic mass conservation constraint, the phasic entropy transport constraint, and the phasic Lin constraint, imposed respectively by the Lagrange multipliers $\xi_k$, $\theta_k$; and $\lambda_{k,i}$. One then considers independent, smooth, and compactly supported perturbations $\delta\mathcal{X}$, for all variables $\mathcal{X} \in \{m_k, \mathbf{u}_k, s_k, \alpha_1, \zeta, \xi_k, \theta_k, \lambda_{k,i}, \mu_{k,i}\}$. The fully Eulerian Stationary Action Principle yields an evolution equation for each variable, except the phasic velocities which satisfy

$$ \mathbf{u}_k = \nabla \xi_k - \theta_k \nabla s_k - \zeta \nabla \alpha_1 - \sum_i \lambda_{k,i} \nabla \mu_{k,i}. \qquad (3.10) $$

Combining the various equations then allows to recover the phasic momentum equations and to derive system (3.6). The details of this derivation are provided in Appendix B.

The decomposition of the phasic velocities given by eq. (3.10) corresponds to a representation in terms of Clebsch potentials, such a decomposition is standard when considering the Stationary Action Principle in an Eulerian framework [36]. The importance of adding Lin's constraints is also highlighted by eq. (3.10). Indeed, assume momentarily that no such constraint was added. Then one obtains (3.10) with $\lambda_{k,i} \equiv 0$ and $\mu_{k,i} \equiv 0$. In particular, in the three-dimensional case, the phasic vorticity then writes

$$ \mathbf{rot}(\mathbf{u}_k) = -\nabla \theta_k \wedge \nabla s_k - \nabla \zeta \wedge \nabla \alpha_1. \qquad (3.11) $$
This is a non-physical and restrictive condition that limits the validity of the equations to a particular subset of flows. For instance, in isentropic flows, it restricts the initial conditions to velocity fields which have their vorticity orthogonal to the volume fraction gradient. Without Lin's constraint in the derivation of the Euler equations, a similar non-physical situation occurs where the vorticity is required to be orthogonal to the gradient of the specific entropy [36].

## 3.2 Analysis of the interfacial quantities

We now discuss the interfacial quantities present in model (3.6) derived in the previous paragraph. The supplementary variable $v$, and the role of the momentum exchange term $f$, will be discussed hereafter.

Recall that model (3.6) relies on three interfacial quantities, namely an interfacial velocity, chosen as the mass-weighted average velocity $u$, an interfacial pressure $p_I$, and an interfacial work term, denoted $(pu)_I$. Although the expression of the interfacial velocity is postulated and imposed as a constraint in the variational principle, the obtained closures for the interfacial pressure (3.4a), and interfacial work (3.7), result from the variational principle itself. They are neither a priori assumptions, nor are they a posteriori closures as in the phase-averaged models. Yet, similar closures may be derived through phase averaging, indeed in [9], a one-dimensional phase-averaged model was derived. The phase-average operator was applied to the entropy equations instead of the energy equations. The interfacial closures compatible with the conservation of the total energy were investigated and it was shown that total energy conservation requires the interfacial velocity and interfacial pressure to satisfy

$$u_1p_1 - u_2p_2 + (u_2 - u_1)p_I + u_I(p_2 - p_1) = 0. \quad (3.12)$$

The authors then proposed $u_I = (u_1 + u_2)/2$, and obtained $p_I = (p_1 + p_2)/2$. However, when inserting $u_I = Y_1u_1 + Y_2u_2$ in (3.12), one obtains $p_I = Y_2p_1 + Y_1p_2$, and thus our model is also compatible with the phase-averaging approach proposed in [9].

The physical interpretation of the interfacial pressure (3.4a) is the following. Each phase is subjected to its own pressure gradient, however, when the phases are out of pressure equilibrium, each phase is subjected to a supplementary force located at the interface and proportional to the exterior pressure. For phase 1, its exterior pressure is $p_2$, whereas for phase 2, its exterior pressure is $p_1$. Averaging the exterior pressures weighted by the concentrations of the phases upon which they act then yields $p_I = Y_1p_2 + Y_2p_1$, i.e. (3.4a), when the concentrations are chosen to be the mass fractions. In the phase averaged model derived in [11], a similar argument was made, but the volume fractions were chosen as concentrations instead, leading to the expression $p_I = \alpha_1p_2 + \alpha_2p_1$. This argument also provides an interpretation of the interfacial work (3.7). If the force induced in phase 1 by the pressure imbalance is proportional to $p_2$, then its associated work is proportional to $p_2u_1$. By symmetry, for phase 2, we obtain $p_1u_2$. Averaging these contributions, in accordance with their respective mass fractions, yields $(pu)_I = Y_1p_2u_1 + Y_2p_1u_2$, which is exactly the expression obtained in (3.7).

Interestingly, our expression for the interfacial pressure (3.4a) is identical to the one obtained in the models derived through the theory of thermodynamically compatible systems, as well as the interfacial pressure (1.1) of the all-topology Baer-Nunziato model analyzed in [14] in the isothermal case. However, as discussed in the introduction, the interfacial pressure (1.1) of the all-topology Baer-Nunziato model strongly depends on the phasic temperatures, while the thermodynamically compatible systems present an artificial heat exchange as well as temperature gradients in their momentum equations. The interfacial work (3.7), provided by the Stationary Action Principle, allows to overcome these issues.

Several other contributions also introduced the interfacial work as an independent quantity in addition to the interfacial velocity and interfacial pressure [11,20,21]. The approaches that have been followed in [11,20], relied on phase-averaging and the interfacial work appeared as one of the unclosed terms, while in [21], a two-pressure two-fluid model has been derived through the
averaging of a stochastic model and its closures depends on the assumptions on the stochastic process. In the stochastic derivation, an interfacial work term, different from the product of the interfacial velocity with the interfacial pressure in the general case, has been determined.

The impact of the interfacial work is highlighted in the phasic internal energy equations, as system (3.6) implies

$$ \alpha_k \rho_k D_t^k e_k = -\alpha_k p_k \nabla \cdot \mathbf{u}_k + p_k (\mathbf{u}_I - \mathbf{u}_k) \cdot \nabla \alpha_k. \quad (3.13a) $$

The first term on the right-hand of corresponds (3.13a) to the classical compression effects induced by a non divergence-free velocity field. The second term is an additional effect induced by a difference between the phasic velocity and the interfacial velocity. When $(\mathbf{u}_I - \mathbf{u}_k) \cdot \nabla \alpha_k > 0$, the kth phase is compressed, while otherwise it is dilatated. For the Baer-Nunziato system, the phasic internal energy equations write

$$ \alpha_k \rho_k D_t^k e_k = -\alpha_k p_k \nabla \cdot \mathbf{u}_k + p_I (\mathbf{u}_I - \mathbf{u}_k) \cdot \nabla \alpha_k. \quad (3.13b) $$

As such, for the Baer-Nunziato models, the rate at which phasic internal energies increase, due to compressions caused by a relative velocity, is given by the interfacial pressure instead of the phasic pressures.

Similarly, for the phasic entropy equations, our model (3.6) results in the transport of the phasic entropies at their own velocity,

$$ \alpha_k \rho_k D_t^k s_k = 0, \quad (3.14a) $$

while for the Baer-Nunziato model, one has

$$ \alpha_k \rho_k D_t^k s_k = \frac{1}{T_k} (p_I - p_k) (\mathbf{u}_I - \mathbf{u}_k) \cdot \nabla \alpha_k. \quad (3.14b) $$

In particular, any all-topology Baer-Nunziato model leads to entropy exchanges between the phases and a conservation law for the total entropy can only be achieved through an interfacial pressure that depends on the phasic temperatures [14,45].

## 3.3 Multi-dimensional effects

We now discuss the forces $f$, defined by (3.8), which appear on the right-hand side of system (3.6), and the new variable $v$. First, in the one-dimensional setting, this force is zero, and equation (3.6h), which governs the evolution of $v$, can be discarded. In the general case, this force may be non-zero, but always satisfies $W \cdot f = 0$, as such it is a force that represents lift. Moreover, because it depends on $\nabla \alpha_1$, this lift force is localized at the interfaces. Its amplitude and direction is partly determined by $v = \nabla \zeta$, which has the dimension of a velocity, and thus $\zeta$ plays the role of a velocity potential. We recall that the evolution equation for $\zeta$, given by (3.3d), writes

$$ \rho D_t \zeta = p_1 - p_2. \quad (3.15) $$

Consequently, for flows which remain at pressure equilibrium, $\zeta$ will be constant along the fictive trajectories defined by the mixture velocity $\mathbf{u}$. If at the initial time it is constant, it remains so and the lift force is zero. As a result, expression (3.8) corresponds to lift induced by pressure non-equilibrium. At this stage, the question of how to initialize $v = \nabla \zeta$, and its physical meaning, remain under investigation.

The appearance of lift forces in two-phase flow models obtained by variational principles has already been noted, see [25] for instance. However, it has been shown in [35], that for a same physical Lagrangian several systems may be derived, which differ only by the expression of their lift forces. In particular, during the variational process, fluid particle trajectories $\mathbf{\Phi}_1$ and $\mathbf{\Phi}_2$, are used to define the variations of the velocities $\mathbf{u}_1$ and $\mathbf{u}_2$. If instead, one operates a change
of variables from ($u_1, u_2$) to ($u, u_1$) or ($u, u_2$), and associates fictive trajectories to these new velocities, then one obtains a new set of equations which differs only by the expression of its lift forces. Consequently, we have used in our derivation only the true fluid particle trajectories associated to $u_1$ and $u_2$.

Lift forces in the phasic momentum equations have also been obtained in the Thermodynamically Compatible models [24–26], their expression however does not require any new variable as the forces write

$$ \pm \rho Y_1 Y_2 [Y_2 \text{rot}(u_1) + Y_1 \text{rot}(u_2)] \wedge W. \qquad (3.16) $$

This expression may be interpreted as the effect of the rotation of one fluid, determined by $\text{rot}(u_k)$, on the other phase $k'$. However, the derivation of these models relies on the change of variables ($u_1, u_2$) to ($u, W$), and the use of a single fictive trajectory associated to the mixture velocity $u$, while the relative velocity $W$ is expressed through the gradient of the Lagrange multiplier which enforces the phasic mass conservation, as shown in [33]. Moreover, in these models, the variational principle imposes $\text{rot}(W) = 0$, i.e. $\text{rot}(u_1) = \text{rot}(u_2)$, and thus the fluids have the same rate of rotation. In light of these elements, and the results shown in [35] concerning the impact of fictive trajectories on lift forces, the validity of either (3.8) or (3.16) remains unclear.

Yet, lift forces do correspond to a physical and observable reality. In the case of disperse flows, a rotating particle in a carrier fluid is subjected to a lift force, called Magnus force, which at low Reynolds numbers writes [46,47],

$$ f_{\text{Magnus}} = \frac{\pi}{8} D^3 \rho_1 \left( \frac{1}{2} \text{rot}(u_1) - \omega_2 \right) \wedge W, \qquad (3.17) $$

where the index 1 represents the carrier fluid and 2 the particle which has diameter $D$, and rotation $\omega_2$. However, in order to obtain such effects, one should need to consider the rotational energy, proportional to $\omega_k^2$, in the Lagrangian. As this energy has been neglected in the two-phase Lagrangian (2.10), neither (3.8) nor (3.8) represent internal degrees corresponding to particle rotations.

In the absence of a clear physical explanation, one possible solution is to neglect the lift forces (3.8) in our model and to discard the supplementary variable $v$ so that the system becomes

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

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

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

$$ \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{I}_d) - p_I \nabla \alpha_1 = 0, \qquad (3.18d) $$

$$ \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{I}_d) + p_I \nabla \alpha_1 = 0, \qquad (3.18e) $$

$$ \partial_t (\alpha_1 \rho_1 E_1) + \nabla \cdot (\alpha_1 \rho_1 H_1) - (p \mathbf{u})_I \cdot \nabla \alpha_1 = 0, \qquad (3.18f) $$

$$ \partial_t (\alpha_2 \rho_2 E_2) + \nabla \cdot (\alpha_2 \rho_2 H_2) + (p \mathbf{u})_I \cdot \nabla \alpha_1 = 0. \qquad (3.18g) $$

We may then conclude that, while our model allows to have both physically coherent closures for the interfacial quantities, and good mathematical properties as will be shown in the Section 4, the modelling of multi-dimensional effects such as lift requires further efforts. We now move on to the next section which concerns the modelling of dissipative source terms.

## 3.4 Modelling source terms

The Hamiltonian principle provides only reversible systems that are without dissipation. Consequently, any dissipative effects must be added a posteriori. Moreover, as the model we have derived allows for full non-equilibrium between the two fluids, similarly to the Baer-Nunziato model, it must be equipped with dissipative source terms in order to provide a meaningful
description of two-phase flows. These source terms must be such that in the absence of exterior perturbations, the system is driven towards equilibrium. Furthermore, they must be in agreement with the Second Principle of Thermodynamics and induce dissipation in the form of entropy production. As the prognostic variables and the total entropy of our model (3.18) are the same as for the Baer-Nunziato model, the same source terms, as determined in [10,13,21,45,48], can be seamlessly considered. Doing so, we obtain

$$ \partial_t \alpha_1 + \mathbf{u} \cdot \nabla \alpha_1 = S_{\text{mec}}, \qquad (3.19\text{a}) $$

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

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

$$ \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{I}_d) - p_I \nabla \alpha_1 = +S_{\text{kin}}, \qquad (3.19\text{d}) $$

$$ \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{I}_d) + p_I \nabla \alpha_1 = -S_{\text{kin}}, \qquad (3.19\text{e}) $$

$$ \partial_t (\alpha_1 \rho_1 E_1) + \nabla \cdot (\alpha_1 \rho_1 H_1) - (p\mathbf{u})_I \cdot \nabla \alpha_1 = +S_{\text{th}} - p_I S_{\text{mec}} + \mathbf{u}_I \cdot S_{\text{kin}}, \qquad (3.19\text{f}) $$

$$ \partial_t (\alpha_2 \rho_2 E_2) + \nabla \cdot (\alpha_2 \rho_2 H_2) + (p\mathbf{u})_I \cdot \nabla \alpha_1 = -S_{\text{th}} + p_I S_{\text{mec}} - \mathbf{u}_I \cdot S_{\text{kin}}, \qquad (3.19\text{g}) $$

where $S_{\text{mec}}$, $S_{\text{kin}}$ and $S_{\text{th}}$ correspond to mechanical, kinematic, and thermal relaxation terms respectively. These source terms can also be added to the system with lift forces (3.6). The total entropy equation then writes

$$ \begin{aligned} \partial_t (\alpha_1 \rho_1 s_1 + \alpha_2 \rho_2 s_2) + \nabla \cdot (\alpha_1 \rho_1 s_1 \mathbf{u}_1 + \alpha_2 \rho_2 s_2 \mathbf{u}_2) = & \left(\frac{1}{T_1} - \frac{1}{T_2}\right) S_{\text{th}} \\ & + \left(\frac{Y_2}{T_1} + \frac{Y_1}{T_2}\right) (\mathbf{u}_2 - \mathbf{u}_1) \cdot \mathbf{S}_{\text{kin}} \quad (3.20) \\ & + \left(\frac{Y_1}{T_1} + \frac{Y_2}{T_2}\right) (p_1 - p_2) S_{\text{mec}}, \end{aligned} $$

and therefore, setting

$$ S_{\text{mec}} = (p_1 - p_2)/\epsilon_p, \qquad (3.21\text{a}) $$

$$ S_{\text{kin}} = (\mathbf{u}_2 - \mathbf{u}_1)/\epsilon_u, \qquad (3.21\text{b}) $$

$$ S_{\text{th}} = (T_2 - T_1)/\epsilon_T, \qquad (3.21\text{c}) $$

then yields a non-negative entropy production for positive relaxation parameters $\epsilon_p$, $\epsilon_u$, and $\epsilon_T$. Several closures for the relaxation parameters have been proposed in the literature, see for instance [21,49–51]. Furthermore, in flow regimes where some of the associated relaxation time-scales are sufficiently small, reduced order models have been derived assuming that the corresponding flow variables remain at equilibrium [48,52–54]. The same model reductions can be applied to (3.19). Moreover, our expression for the interfacial work (3.7) satisfies

$$ (p\mathbf{u})_I = p_I \mathbf{u}_I - Y_1 Y_2 (\mathbf{u}_1 - \mathbf{u}_2) (p_1 - p_2), \qquad (3.22) $$

and therefore, when considering the asymptotic expansions near velocity or pressure equilibrium, the interfacial work represents a first order perturbation of a Baer-Nunziato model, and thus the reduced order models obtained at zeroth order are identical. This is a reassuring result for two reasons. First, the closure problem for interfacial quantities is related to non-equilibrium flows only. Second, several of these reduced order models recover some well-known two-phase flow features, such as the so-called Kapila model [48], which is derived assuming velocity and pressure equilibrium, and recovers the non-monotonic mixture sound speed as determined by Wood [55].

Unfortunately, unless complete equilibrium is assumed, these reduced order models are known to be either non-hyperbolic [54], or to have non-conservative products which are ill-defined for shocks [48]. Therefore, the use of fully out of equilibrium models, such as the one
we have derived, combined with relaxation source terms and numerical methods capable of handling their inherent stiffness, such as the one proposed in [51], is preferable to obtain convergent numerical simulations.

The next section is dedicated to the mathematical analysis of the model, in particular we exhibit the necessary properties for its well-posedness and show that it allows for the proper treatment of weak solutions.

# 4 Mathematical analysis of the model: the one-dimensional setting

In this section, the mathematical properties of the model are analysed. We refer to textbooks on systems of conservation laws, such as [56], for definitions of hyperbolicity, Riemann invariants, and the various properties that will be considered. Our analysis will be carried-out in the one-dimensional setting and without source terms. In this case, the lift forces are zero and the equation on $v$ is discarded, as in (3.18). As the system is close to the Baer-Nunziato model, many properties verified for the Baer-Nunziato model are also verified for our model. Moreover, the proofs provided are adaptations of existing proofs for the Baer-Nunziato model found in [14,16]. Some properties concerning the multi-dimensional case are provided at the end of this section.

## 4.1 Hyperbolicity and characteristic fields

The development of two-pressure models is mainly motivated by the ill-posedness of the single-pressure models [57]. Single-pressure models exhibit complex characteristics and are therefore linearly unstable. Consequently, we consider real characteristics to be a minimal requirement for the modelling of two-phase flows. However, several contributions have suggested that instabilities due to non-hyperbolicity may be controlled [6,7]. Two-pressure models often have real characteristics [9,14]. For the model we have derived, this is specified by the following proposition.

**Proposition 4.1 (Hyperbolicity).** The system admits real eigenvalues which are $u, u_k$ and $u_k \pm c_k$, with $k = 1, 2$. The corresponding eigenvectors span $\mathbb{R}^7$, and the system is hyperbolic, under the non-resonance condition $u \neq u_k \pm c_k$, for $k = 1, 2$.

*Proof.* Through a change of variables, the system is expressed in terms of the variable $U = (\alpha_1, u_1, p_1, s_1, u_2, p_2, s_2)^T$, and writes in the following quasilinear form,

$$ \partial_t U + \mathbf{M}(U)\partial_x U = 0, \quad (4.1) $$

with

$$ \mathbf{M}(U) = \begin{bmatrix} u & \mathbf{0}_{1 \times 3} & \mathbf{0}_{1 \times 3} \\ \mathbf{z}_{\alpha,1} & \mathbf{M}_1 & \mathbf{0}_{3 \times 3} \\ \mathbf{z}_{\alpha,2} & \mathbf{0}_{3 \times 3} & \mathbf{M}_2 \end{bmatrix}, \quad (4.2) $$

and

$$ \mathbf{M}_k = \begin{bmatrix} u_k & 1/\rho_k & 0 \\ \rho_k c_k^2 & u_k & 0 \\ 0 & 0 & u_k \end{bmatrix}, \quad \mathbf{z}_{\alpha,k} = \begin{bmatrix} \frac{p_1-p_2}{\rho_k} \\ \frac{\rho_k c_k^2}{\alpha_k} (u-u_k) \\ 0 \end{bmatrix}. \quad (4.3) $$

The block-triangular structure of $\mathbf{M}(U)$ allows for an immediate computation of the eigenvalues. When the non-resonance condition is satisfied, the matrix $\mathbf{M}(U)$ is diagonalizable and a matrix
of right eigenvectors $\mathbf{R}(U)$ is given by

$$ \mathbf{R}(U) = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ \frac{-Y_2 W}{(u_1-u)^2 - c_1^2} \left( \frac{p_1-p_2}{\rho} - \frac{c_1^2}{\alpha_1} \right) & 1 & 1 & 0 & 0 & 0 & 0 \\ \frac{\rho_1 c_1^2}{(u_1-u)^2 - c_1^2} \left( \frac{p_1-p_2}{\rho} - \frac{Y_2^2 W^2}{\alpha_1} \right) & \rho_1 c_1 & -\rho_1 c_1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ \frac{Y_1 W}{(u_2-u)^2 - c_2^2} \left( \frac{p_1-p_2}{\rho} + \frac{c_2^2}{\alpha_2} \right) & 0 & 0 & 0 & 1 & 1 & 0 \\ \frac{\rho_2 c_2^2}{(u_2-u)^2 - c_2^2} \left( \frac{p_1-p_2}{\rho} + \frac{Y_1^2 W^2}{\alpha_2} \right) & 0 & 0 & 0 & \rho_2 c_2 & -\rho_2 c_2 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{bmatrix}. \quad (4.4) $$

The columns of $\mathbf{R}(U)$ correspond to eigenvectors $u, u_1+c_1, u_1-c_1, u_1, u_2+c_2, u_2-c_2$ and $u_2$ respectively. □

**Remark 4.1.** The non-resonance condition also applies to the Baer-Nunziato models [14] and the two-pressure models derived in [9].

**Remark 4.2.** At resonance, two eigenvalues coincide and an eigenvector may be missing. However, as determined in [9], a complete set of eigenvectors may be available under additional conditions. Assume resonance for the kth phase, i.e. $c_k^2 = (u - u_k)^2$, if $p_I = p_k - \rho_k c_k^2$, then the system remains hyperbolic. This also holds when both phases are in resonance.

We now investigate the nature of the different characteristic fields. To do so, we introduce, for each phase $k$, its so-called fundamental derivative $G_k$ [58], which is the thermodynamic function defined by

$$ G_k = -\frac{\tau_k \partial^3 e_k / \partial \tau_k^3 |_{s_k}}{2 \partial^2 e_k / \partial \tau_k^2 |_{s_k}}, \quad (4.5) $$

where $\tau_k = 1/\rho_k$, is the phase's specific volume.

The fundamental derivative $G_k$ plays an important role in the study of the Riemann problem and the more restrictive condition $G_k > 0$, or even $G_k > 1$, is usually implicitly assumed. For a single fluid, the condition $G > 0$ corresponds to the convexity of the isentropes in the $p-\tau$ plane, with $\tau = 1/\rho$, the specific volume. In particular, it ensures that shocks are compressive, while for $G < 0$, the system admits rarefactions shocks instead. The stronger condition, $G > 1$, expresses the convexity of the isentropes in the $p-\rho$ plane and ensures that the sound velocity increases along the Hugoniot curves [58,59].

We then have the following classification of the characteristic fields.

**Proposition 4.2 (Characteristic fields).** When the system is out of resonance, the characteristic fields associated to the convective eigenvalues $u$ and $u_k$, for $k = 1, 2$, are linearly degenerate. Moreover, for $k = 1, 2$, if the equation of state of the kth phase is such that $G_k \neq 0$, then the characteristic fields associated to the acoustic eigenvalues $u_k \pm c_k$, are genuinely non-linear.

*Proof.* The linearly degenerate character of the eigenvalues $u_1$ and $u_2$ is immediate owing to the expression of their eigenvectors given by (4.4). For $k = 1, 2$, let $\lambda_k(U) = u_k \pm c_k$ and $r_k(U)$ be the associated eigenvector as given by (4.4). We then have that

$$ \nabla_U \lambda_k \cdot r_k = 1 + \rho_k c_k \frac{\partial c_k}{\partial p_k|_{s_k}} = 1 + \frac{1}{2} \rho_k \frac{\frac{\partial^2 p_k}{\partial \rho_k^2}|_{s_k}}{\frac{\partial p_k}{\partial \rho_k}|_{s_k}} = G_k, \quad (4.6) $$

which shows the genuinely non-linear character of the acoustic waves. As the classification of the characteristic fields does not depend on the choice of variables, to show the linear degeneracy of the interfacial characteristic field, we use the set of variables $V = (\alpha_1, u_1, p_1, m_1, u_2, p_2, m_2)$
with $m_k = \alpha_k \rho_k$. With these variables, the eigenvector $r$ associated to the eigenvalue $u$ is given by the first column of (4.4) but where its two zero entries are respectively replaced by $\frac{\rho Y_1}{(u_1-u)^2-c_1^2}$, $((p_1-p_2)/\rho - c_1^2/\alpha_1)$ and $\frac{\rho Y_2}{(u_2-u)^2-c_2^2}$, $((p_1-p_2)/\rho + c_2^2/\alpha_2)$. Moreover, since $\nabla_V u = (0, Y_2 W/\rho, Y_1, 0, -Y_1 W/\rho, Y_2, 0)^\top$, a direct computation yields $\nabla_V u \cdot r = 0$. $\square$

## 4.2 Riemann invariants and jump conditions

Solutions to the Riemann Problem consist of constant states connected through either contact discontinuities associated to linearly degenerate waves, or shock or rarefaction waves associated to genuinely non-linear waves. States that can be connected to each other through contact discontinuities or rarefactions waves can be parametrized by the Riemann invariants of the corresponding wave. Through such waves, the Riemann invariants remain constant.

**Proposition 4.3 (Riemann invariants).** *Outside of resonance, each characteristic field admits a set of 6 Riemann invariants with linearly independent gradients. For each characteristic field, such a set of Riemann invariants is given by*

$$u: \quad \{u, \rho Y_1 Y_2 W, \alpha_1 p_1 + \alpha_2 p_2 + \rho Y_1 Y_2 W^2, h_1 - h_2 + \frac{1}{2}(Y_2 - Y_1)W^2, s_1, s_2\}, \quad (4.7a)$$

$$u_1: \quad \{\alpha_1, u_1, p_1, u_2, p_2, \rho_2\}, \quad (4.7b)$$

$$u_1 + c_1: \quad \{\alpha_1, s_1, u_1 - f_1(s_1, p_1), u_2, p_2, \rho_2\}, \quad (4.7c)$$

$$u_1 - c_1: \quad \{\alpha_1, s_1, u_1 + f_1(s_1, p_1), u_2, p_2, \rho_2\}, \quad (4.7d)$$

$$u_2: \quad \{\alpha_1, u_2, p_2, u_1, p_1, \rho_1\}, \quad (4.7e)$$

$$u_2 + c_2: \quad \{\alpha_1, s_2, u_2 - f_2(s_2, p_2), u_1, p_1, \rho_1\}, \quad (4.7f)$$

$$u_2 - c_2: \quad \{\alpha_1, s_2, u_2 + f_2(s_2, p_2), u_1, p_1, \rho_1\}, \quad (4.7g)$$

where $f_k(s_k, p_k) = \int^{p_k} [\rho_k(s_k, p)c_k(s_k, p)]^{-1} dp$, for $k = 1, 2$.

*Proof.* The Riemann invariants are obtained by direct computation using the expression of the eigenvectors (4.4). For some of the Riemann invariants associated to the interface wave (4.7a), it is more practical to use the variables $V$, and the associated eigenvector, introduced in the proof of proposition 4.2, in combination with the following thermodynamic identity

$$c_k^2 \frac{\partial h_k}{\partial p_k} + \frac{\partial h_k}{\partial \rho_k} = \frac{c_k^2}{\rho_k}. \quad (4.8)$$

**Remark 4.3.** *The interfacial pressure $p_I$, given by eq. (3.4a), is not continuous through the interfacial wave. Instead it is a total pressure $\alpha_1 p_1 + \alpha_2 p_2 + \rho Y_1 Y_2 W^2$ which is continuous. This expression was chosen as interfacial pressure for a Baer-Nunziato model in [13], but leads to a non-conservative equation for the total entropy equation and thus the model lacks an admissibility criterion for shocks.*

The Riemann invariants are identical to the ones of the Baer-Nunziato system, see [14], except for those associated to the interfacial wave $u$.

In addition, similarly to the Baer-Nunziato models, the two-phases are coupled only through the non-conservative products which involve $\partial_x \alpha_1$. Owing to its linear degeneracy, the non-conservative products through the interfacial wave $u$ are entirely defined by the Riemann invariants (4.7a). Moreover, the volume fraction remains constant on both sides of the interfacial wave, and thus the non-conservative products are inactive through shocks and the phases are decoupled [14]. In particular, if we denote by $[\phi]$ the jump of a variable $\phi$ through a discontinuity, we obtain the following jump conditions for the genuinely non-linear waves.
**Proposition 4.4 (Jump conditions).** For the genuinely non-linear waves, two states connected through a shock, propagating at speed $\sigma$, satisfy the following jump conditions,

$$ \llbracket \alpha_1 \rrbracket = 0, \qquad (4.9a) $$

$$ \llbracket \alpha_k \rho_k (u_k - \sigma) \rrbracket = 0, \qquad (4.9b) $$

$$ \llbracket \alpha_k \rho_k u_k (u_k - \sigma) \rrbracket + \llbracket \alpha_k p_k \rrbracket = 0, \qquad (4.9c) $$

$$ \llbracket \alpha_k \rho_k E_k (u_k - \sigma) \rrbracket + \llbracket \alpha_k p_k u_k \rrbracket = 0, \qquad (4.9d) $$

for $k = 1, 2$.

*Proof.* Direct application of the Rankine-Hugoniot jump conditions when $\partial_x \alpha_1 = 0$. $\square$

In Section 3, we have imposed the transport of specific entropies for each phase as a constraint in the variational derivation of the model. Consequently, the model admits a supplementary conservation law for the physical entropy which allows to select physically admissible shocks. Note that the availability of an entropy conservation law for the Baer-Nunziato models is not guaranteed and depends on the closures for the interfacial quantities [14].

A shock is said to be admissible if it satisfies the following entropy inequality,

$$ \partial_t (\alpha_1 \rho_1 s_1 + \alpha_2 \rho_2 s_2) + \partial_x (\alpha_1 \rho_1 s_1 u_1 + \alpha_2 \rho_2 s_2 u_2) \geq 0, \qquad (4.10) $$

and thus the mathematical entropy is given by the opposite of the physical entropy. Note that when written in terms of mixture velocity $u$ and relative velocity $W$, eq. (4.10) corresponds to (1.2).

**Remark 4.4.** The mathematical entropy is convex, but not strictly convex as shown in [16].

## 4.3 Symmetrization

The previous results concern weak solutions, and in particular, we have shown that shocks are uniquely defined despite the non-conservative nature of the system. Here, we tackle the question of its well-posedness for smooth initial data, in the sense of Sobolev regularity. Friedrichs introduced an important class of symmetric hyperbolic systems [60], for which local-in-time well-posedness results from Kato's theorem [61]. These results then extend to symmetrizable systems.

In [16], the Baer-Nunziato models have been shown to be symmetrizable when out of resonance, independently of the closures for the interfacial quantities. Although the derived model is not a Baer-Nunziato model, it is, in some sense, an extension of this class of models and thus we have the following proposition.

**Proposition 4.5 (Symmetrizability).** The system is symmetrizable when the non-resonance condition is satisfied.

*Proof.* We follow and adapt the proof of the symmetrizability of the Baer-Nunziato system proposed in [16]. We express the system in the same quasilinear form (4.1) as in the proof of proposition 4.1 and exhibit a symmetric positive definite matrix $P(U)$ such that $P(U)M(U)$ is symmetric. A symmetric form of the system is then given by

$$ P(U) \partial_t U + P(U) M(U) \partial_x U = 0. \qquad (4.11) $$

Setting $P_k = \text{diag}(\rho_k c_k^2, 1, 1)$, we consider

$$ P = \begin{bmatrix} \theta_{\alpha} & y_{1,\alpha}^{\top} & y_{2,\alpha}^{\top} \\ y_{1,\alpha} & P_1 & 0_{3 \times 3} \\ y_{2,\alpha} & 0_{3 \times 3} & P_2 \end{bmatrix}, \qquad (4.12) $$
where $\theta_{\alpha} \in \mathbb{R}$ and $\mathbf{y}_{1,\alpha}, \mathbf{y}_{2,\alpha} \in \mathbb{R}^3$ are to be determined. The requirement $\mathbf{P}(U)\mathbf{M}(U)$ symmetric is equivalent to

$$ (\mathbf{M}_k^\top - u\mathbf{I}_3) \mathbf{y}_{k,\alpha} = \mathbf{P}_k \mathbf{z}_{k,\alpha}, \quad k = 1, 2. \qquad (4.13) $$

This is a 3 × 3 linear system whose last line yields a trivial equation and allows to set $\mathbf{y}_{k,\alpha}^{(3)} = 0$, the remaining 2 × 2 system writes

$$ \begin{bmatrix} u_k - u & \rho_k c_k^2 \\ 1/\rho_k & u_k - u \end{bmatrix} \begin{bmatrix} \mathbf{y}_{k,\alpha}^{(1)} \\ \mathbf{y}_{k,\alpha}^{(2)} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix} \mathbf{M}_k \mathbf{z}_{k,\alpha}, \quad k = 1, 2. \qquad (4.14) $$

When the non-resonance condition is satisfied, this linear system is invertible and yields a closed expression for $\mathbf{y}_{k,\alpha}$, for $k = 1, 2$, such that $\mathbf{P}(U)\mathbf{M}(U)$ is symmetric. It remains to determine $\theta_{\alpha}$ such that $\mathbf{P}(U)$ is positive definite. For any $\mathbf{X} = (a, \mathbf{b}_1^{\top}, \mathbf{b}_2^{\top}) \in \mathbb{R}^7$, we have that

$$ \mathbf{X}^{\top} \mathbf{P}(U) \mathbf{X} = a^2 \theta_{\alpha} + 2a (\mathbf{y}_{1,\alpha}^{\top} \mathbf{b}_1 + \mathbf{y}_{2,\alpha}^{\top} \mathbf{b}_2) + (\mathbf{b}_1^{\top} \mathbf{P}_1 \mathbf{b}_1 + \mathbf{b}_2^{\top} \mathbf{P}_2 \mathbf{b}_2), \qquad (4.15) $$

This is a second-order polynomial with respect to $a$ and its discriminant $\Delta$ writes

$$ \Delta = 4 \left[ (\mathbf{y}_{1,\alpha}^{\top} \mathbf{b}_1 + \mathbf{y}_{2,\alpha}^{\top} \mathbf{b}_2)^2 - \theta_{\alpha} (\mathbf{b}_1^{\top} \mathbf{P}_1 \mathbf{b}_1 + \mathbf{b}_2^{\top} \mathbf{P}_2 \mathbf{b}_2) \right]. \qquad (4.16) $$

As $\mathbf{P}_k$ is symmetric positive definite, its square root, $\mathbf{P}_k^{1/2}$, is also symmetric positive definite and we may set $\bar{\mathbf{b}}_k = \mathbf{P}_k^{1/2} \mathbf{b}_k$, for $k = 1, 2$, such that

$$ \frac{\Delta}{4} = (\mathbf{y}_{1,\alpha}^{\top} \mathbf{P}_1^{-1/2} \bar{\mathbf{b}}_1 + \mathbf{y}_{2,\alpha}^{\top} \mathbf{P}_2^{-1/2} \bar{\mathbf{b}}_2)^2 - \theta_{\alpha} (|\bar{\mathbf{b}}_1|^2 + |\bar{\mathbf{b}}_2|^2). \qquad (4.17) $$

Applying the Cauchy-Schwarz inequality in $\mathbb{R}^6$ to $(\bar{\mathbf{b}}_1^{\top}, \bar{\mathbf{b}}_2^{\top})$ and $((\mathbf{P}_1^{-1/2} \mathbf{y}_{1,\alpha})^{\top}, (\mathbf{P}_2^{-1/2} \mathbf{y}_{2,\alpha})^{\top})$ then yields

$$ \frac{\Delta}{4} \leq \left( |\mathbf{P}_1^{-1/2} \mathbf{y}_{1,\alpha}|^2 + |\mathbf{P}_2^{-1/2} \mathbf{y}_{2,\alpha}|^2 - \theta_{\alpha} \right) (|\bar{\mathbf{b}}_1|^2 + |\bar{\mathbf{b}}_2|^2). \qquad (4.18) $$

Choosing $\theta_{\alpha} > |\mathbf{P}_1^{-1/2} \mathbf{y}_{1,\alpha}|^2 + |\mathbf{P}_2^{-1/2} \mathbf{y}_{2,\alpha}|^2$ results in a non-positive discriminant so that $\mathbf{X}^{\top} \mathbf{P}(U) \mathbf{X} \ge 0$. Moreover, if $\mathbf{X}^{\top} \mathbf{P}(U) \mathbf{X} = 0$, then necessarily $\Delta = 0$ and so by (4.18), $\bar{\mathbf{b}}_k = 0$. It is then immediate that $\mathbf{b}_k = 0$, and by (4.15) it follows that $a = 0$, so that $\mathbf{X} = 0$. As a result, we may choose $\theta_{\alpha}$ such that $\mathbf{P}(U)$ is symmetric positive definite and symmetrizes the quasilinear form (4.1) of the system. $\square$

**Remark 4.5.** For conservative systems, the existence of a supplementary conservation law for a strictly convex entropy implies symmetrizability of the system [62,63]. The link between the existence of a supplementary conservation law and symmetrizability for non-conservative systems has been investigated in [64]. In particular, the proposed framework of [64] requires the strict convexity of the supplementary potential only for a subset of the system's variables and has been applied to the Baer-Nunziato system.

As a result of its symmetrizability, for initial conditions far from resonance, the Cauchy problem admits a local-in-time smooth solution [61]. The smooth solution exists until it either becomes resonant or has derivatives which become unbounded in $L^{\infty}$.
## 4.4 Some comments on the multi-dimensional case

Owing to the model's rotational invariance, its hyperbolicity in the multi-dimensional case may studied only in the *x*-direction. If the lift forces are neglected and one considers system (3.18), then the model remains hyperbolic. It has the same eigenvalues as in the one-dimensional case but the multiplicity of the convective waves $u_{1,x}$ and $u_{2,x}$ is increased by the number of tangential components. Moreover, the nature of the fields remains identical and thus the previous discussion concerning shocks and non-conservative products remains valid.

If one chooses to keep the lift forces and considers system (3.6), then the model is only weakly hyperbolic. Again, the model has the same eigenvalues as in the one-dimensional case, and the multiplicity of the convective waves $u_{1,x}$ and $u_{2,x}$ is also increased by the number of tangential components. However, owing to the additional equation (3.6h) for the supplementary variable $v$, the multiplicity of the interfacial eigenvalue $u_x$ is also increased but one eigenvector is missing and we do not have a full basis of eigenvectors.

This is similar to the Thermodynamically Compatible two-phase models [25,26], which are also weakly hyperbolic in the multi-dimensional case. In fact, weak hyperbolicity is recurrent for systems which have an involution constraint such as $\text{rot}(v) = 0$, see for instance [26,40]. Several techniques to restore hyperbolicity have been developed, such as curl-cleaning [65], and Godunov-Powell symmetrization [66,67] but have not been investigated in the present work.

# 5 Conclusion

In this paper, the generic Stationary Action Principle framework presented in [39] has been extended to the case of two families of trajectories. The new framework allows for the derivation of the equations of motion in two-velocity models for two-phase flows and guarantees the conservation of total momentum and total energy.

The newly developed framework is then used for the derivation of a novel *all-topology* two-fluid model. The model is reminiscent of the class of Baer-Nunziato models that it extends through the introduction of a new interfacial quantity, namely the interfacial work. However, contrary to the Baer-Nunziato models, for which many different closures have been considered, our model is fully closed as the sole assumption on the interfacial velocity, in the variational principle, suffices to close the remaining interfacial quantities. The interfacial velocity is given by the mixture velocity, i.e. the mass-weighted average velocity, allowing it to be well-defined and have a meaningful expression regardless of flow conditions. Moreover, the resulting expressions for the interfacial pressure and interfacial work are physically sound and allow to clearly distinguish the various mechanical effects from the thermodynamic ones.

In the one-dimensional setting, the model is shown to be hyperbolic under a classical non-resonance condition and is symmetrizable, allowing for local-in-time well-posedness for smooth solutions. Moreover, despite its non-conservative nature, its non-conservative products are uniquely defined through discontinuities and an entropy conservation law provides a selection criteria for physically admissible shocks.

In the multi-dimensional setting, lift forces are present in the phasic momentum equations and are driven by a new potential velocity which has its own evolution equation and the resulting model is weakly hyperbolic. However, in the absence of rotational energies in the Lagrangian, a clear interpretation of the obtained lift forces is not yet available and requires further investigation. Consequently, we may choose to neglect these forces at the present time, and in doing so, we obtain a model which inherits the mathematical properties proved in the one-dimensional case.

It is, to the best of our knowledge, the first *all-topology* multi-fluid model that has both a complete physical consistency and the required mathematical structure to properly describe shock-type solutions. The contribution of this paper is the development of a model, but this model lends itself well to simulations even though this is not within the scope of the paper.
Readers can find some simple numerical simulations at the address indicated<fnref n="1" />, along with the codes to reproduce them. Future works will concern the development of accurate approximate Riemann solvers for the model and further investigations in the modelling of lift forces.

**Acknowledgements.** The authors would like to thank Jean-Marc Hérard for fruitful discussions on related topics. This work has been supported by the Agence Innovation Défense (AID), the CIEDS project OPEN-NUM-DEF (PI. L. Gouarin, M. Massot, T. Pichard), a PhD Grant from ONERA (W. Haegeman), as well as the HPC@Maths Initiative (PI. L. Gouarin and M. Massot) of the Fondation École polytechnique, which we would like to acknowledge.

# A Derivation of the generic system through Stationary Action principle

## A.1 Proof of Theorem 2.1

We provide here a proof of the generic equations (2.5) obtained by the Stationary Action Principle (2.2), for the generic Lagrangian (2.1). The variation of the action $\mathcal{A}$ is given by

$$ \delta \mathcal{A} = \sum_{k=1}^{2} \iint_{\Omega_t} \left\{ \frac{\partial \mathcal{L}_k}{\partial m_k} \delta m_k + \frac{\partial \mathcal{L}_k}{\partial \mathbf{u}_k} \delta \mathbf{u}_k + \sum_{i=1}^{n_k} \frac{\partial \mathcal{L}_k}{\partial b_{k,i}} \delta b_{k,i} + \sum_{i=1}^{n_{k'}} \frac{\partial \mathcal{L}_k}{\partial b_{k',i}} \delta b_{k',i} \right\} d\mathbf{x}dt \\ + \sum_{k=1}^{2} \sum_{j=1}^{m} \iint_{\Omega_t} \left\{ \frac{\partial \mathcal{L}_k}{\partial \xi_j} \delta \xi_j + \frac{\partial \mathcal{L}_k}{\partial (\nabla \xi_j)} \delta (\nabla \xi_j) + \frac{\partial \mathcal{L}_k}{\partial (D_t^k \xi_j)} \delta (D_t^k \xi_j) \right\} d\mathbf{x}dt. \quad (A.1) $$

Here and in the sequel, $k' = 3 - k$ and is the opposite phase of phase $k$. For the sake of convenience on the notation, the differential element $d\mathbf{x}dt$ will be dropped in the remaining of the proof. The perturbations appearing in (A.1) do not form free family as they are related to the independent perturbations $\delta\mathbf{\Phi}_1$, $\delta\mathbf{\Phi}_2$ and $\delta\xi_j$ through eqs. (2.4). Using these expressions, the various terms in the previous equation write,

$$ \iint_{\Omega_t} \frac{\partial \mathcal{L}_k}{\partial m_k} \delta m_k = - \iint_{\Omega_t} \frac{\partial \mathcal{L}_k}{\partial m_k} \nabla \cdot (m_k \delta \mathbf{\Phi}_k) = \iint_{\Omega_t} \delta \mathbf{\Phi}_k^\top m_k \nabla \left( \frac{\partial \mathcal{L}_k}{\partial m_k} \right), \quad (A.2a) $$

$$ \iint_{\Omega_t} \frac{\partial \mathcal{L}_k}{\partial b_{k,i}} \delta b_{k,i} = - \iint_{\Omega_t} \delta \mathbf{\Phi}_k^\top \frac{\partial \mathcal{L}_k}{\partial b_{k,i}} \nabla b_{k,i}, \quad (A.2b) $$

$$ \iint_{\Omega_t} \frac{\partial \mathcal{L}_k}{\partial b_{k',i}} \delta b_{k',i} = - \iint_{\Omega_t} \delta \mathbf{\Phi}_{k'}^\top \frac{\partial \mathcal{L}_k}{\partial b_{k',i}} \nabla b_{k',i}, \quad (A.2c) $$

and

$$ \begin{aligned} \iint_{\Omega_t} \frac{\partial \mathcal{L}_k}{\partial \mathbf{u}_k} \delta \mathbf{u}_k &= \iint_{\Omega_t} \frac{\partial \mathcal{L}_k}{\partial \mathbf{u}_k} \left[ D_t^k (\delta \mathbf{\Phi}_k) - (\partial_x \mathbf{u}_k) \delta \mathbf{\Phi}_k \right] \\ &= - \iint_{\Omega_t} \delta \mathbf{\Phi}_k^\top \left[ \partial_t \mathbf{K}_k + \nabla \cdot (\mathbf{K}_k \otimes \mathbf{u}_k) + (\nabla \mathbf{u}_k) \mathbf{K}_k \right]. \end{aligned} \quad (A.2d) $$

Recall that the perturbations are assumed to vanish at the boundary and therefore no boundary integral contribution arises. The independent variables contribute through the following terms

$$ \iint_{\Omega_t} \frac{\partial \mathcal{L}_k}{\partial (\nabla \xi_j)} \delta (\nabla \xi_j) = \iint_{\Omega_t} \frac{\partial \mathcal{L}_k}{\partial (\nabla \xi_j)} \nabla (\delta \xi_j) = - \iint_{\Omega_t} \delta \xi_j \nabla \cdot (\mathbf{D}_{k,j}), \quad (A.3) $$

<footnote n="1">¹https://hpc-maths.github.io/2025_09_two_fluid_all_topology/</footnote>
and

$$
\begin{align}
\iint_{\Omega_t} \frac{\partial \mathcal{L}_k}{\partial (\mathrm{D}_t^k \xi_j)} \delta(\mathrm{D}_t^k \xi_j) &= \iint_{\Omega_t} \mathcal{M}_{k,j} \mathrm{D}_t^k (\delta \xi_j) \nonumber \\
&\quad + \iint_{\Omega_t} \mathcal{M}_{k,j} \nabla \xi_j^\top [\mathrm{D}_t^k (\delta \mathbf{\Phi}_k) - (\partial_x \mathbf{u}_k) \delta \mathbf{\Phi}_k] \nonumber \\
&= - \iint_{\Omega_t} \delta \xi_j [\partial_t \mathcal{M}_{k,j} + \nabla \cdot (\mathcal{M}_{k,j} \mathbf{u}_k)] \tag{A.4} \\
&\quad - \iint_{\Omega_t} \delta \mathbf{\Phi}_k^\top \Biggl[ \begin{aligned}[t]
    &\partial_t (\mathcal{M}_{k,j} \nabla \xi_j) + \nabla \cdot (\mathcal{M}_{k,j} \nabla \xi_j \otimes \mathbf{u}_k) \\
    &\quad + (\nabla \mathbf{u}_k) \mathcal{M}_{k,j} \nabla \xi_j \end{aligned} \Biggr]. \nonumber
\end{align}
$$

Combining these expressions, we get

$$
\begin{equation}
\begin{split}
\delta \mathcal{A} = - \sum_{k=1}^{2} \iint_{\Omega_t} \delta \mathbf{\Phi}_k^\top & \left\{ \begin{aligned}[t]
    & \partial_t \mathbf{K}_k + \nabla \cdot (\mathbf{K}_k \otimes \mathbf{u}_k) + (\nabla \mathbf{u}_k) \mathbf{K}_k + \sum_{i=1}^{n_k} \left[ \frac{\partial \mathcal{L}_k}{\partial b_{k,i}} + \frac{\partial \mathcal{L}_{k'}}{\partial b_{k,i}} \right] \nabla b_{k,i} \\
    & - m_k \nabla \left( \frac{\partial \mathcal{L}_k}{\partial m_k} \right) + \sum_{j=1}^{m} \left[ \partial_t (\mathcal{M}_{k,j} \nabla \xi_j) + \nabla \cdot (\mathcal{M}_{k,j} \nabla \xi_j \otimes \mathbf{u}_k) + (\nabla \mathbf{u}_k) \mathcal{M}_{k,j} \nabla \xi_j \right]
    \end{aligned} \right\} \\
& - \sum_{j=1}^{m} \iint_{\Omega_t} \delta \xi_j \left\{ \begin{aligned}[t]
    & \partial_t (\mathcal{M}_{1,j} + \mathcal{M}_{2,j}) + \nabla \cdot (\mathcal{M}_{1,j} \mathbf{u}_1 + \mathcal{M}_{2,j} \mathbf{u}_2 + \mathbf{D}_{1,j} + \mathbf{D}_{2,j}) \\
    & - \frac{\partial \mathcal{L}_1}{\partial \xi_j} - \frac{\partial \mathcal{L}_2}{\partial \xi_j}
    \end{aligned} \right\}
\end{split}
\tag{A.5}
\end{equation}
$$

The critical points of the Action then satisfy $\delta\mathcal{A} = 0$, for any compactly supported perturbations $\delta\mathbf{\Phi}_1$, $\delta\mathbf{\Phi}_2$ and $\delta\xi_j$. In particular, the perturbations $\delta\xi_j$ yield the generic equations (2.5c) for $j = 1, \dots, m$. The perturbations $\delta\mathbf{\Phi}_k$, $k = 1, 2$, yield

$$
\begin{equation}
\begin{split}
\partial_t \mathbf{K}_k + \nabla (\mathbf{K}_k \otimes \mathbf{u}_k) + (\nabla \mathbf{u}_k) \mathbf{K}_k - m_k \nabla \left( \frac{\partial \mathcal{L}_k}{\partial m_k} \right) &+ \sum_{i=1}^{n_k} \left[ \frac{\partial \mathcal{L}_k}{\partial b_{k,i}} + \frac{\partial \mathcal{L}_{k'}}{\partial b_{k,i}} \right] \nabla b_{k,i} \\
&+ \sum_{j=1}^{m} \left[ \partial_t (\mathcal{M}_{k,j} \nabla \xi_j) + \nabla \cdot (\mathcal{M}_{k,j} \nabla \xi_j \otimes \mathbf{u}_k) + (\nabla \mathbf{u}_k) \mathcal{M}_{k,j} \nabla \xi_j \right] = 0. \quad (A.6)
\end{split}
\end{equation}
$$

The sum indexed by $j$ on the second line of the previous relation simplifies into

$$
\sum_{j=1}^{m} \left\{ \left[ \partial_t \mathcal{M}_{k,j} + \nabla \cdot (\mathcal{M}_{k,j} \mathbf{u}_k) \right] \nabla \xi_j + \mathcal{M}_{k,j} \nabla \left( \mathrm{D}_t^k \xi_j \right) \right\}. \quad (A.7)
$$

Moreover, we have that

$$
\begin{equation}
\begin{split}
\nabla \mathcal{L}_k = \frac{\partial \mathcal{L}_k}{\partial m_k} \nabla m_k &+ (\nabla \mathbf{u}_k) \mathbf{K}_k + \sum_{i=1}^{n_k} \frac{\partial \mathcal{L}_k}{\partial b_{k,i}} \nabla b_{k,i} + \sum_{i=1}^{n_{k'}} \frac{\partial \mathcal{L}_k}{\partial b_{k',i}} \nabla b_{k',i} \\
&+ \sum_{j=1}^{m} \left\{ \frac{\partial \mathcal{L}_k}{\partial \xi_j} \nabla \xi_j + (\nabla \nabla \xi_j) \mathbf{D}_{k,j} + \mathcal{M}_{k,j} \nabla (\mathrm{D}_t^k \xi_j) \right\}.
\end{split}
\tag{A.8}
\end{equation}
$$
Plugging these expressions into (A.6), we obtain

$$ \partial_t \mathbf{K}_k + \nabla (\mathbf{K}_k \otimes \mathbf{u}_k) + \nabla \left( \mathcal{L}_k - m_k \frac{\partial \mathcal{L}_k}{\partial m_k} \right) - \sum_{i=1}^{n_{k'}} \frac{\partial \mathcal{L}_k}{\partial b_{k',i}} \nabla b_{k',i} + \sum_{i=1}^{n_k} \frac{\partial \mathcal{L}_{k'}}{\partial b_{k,i}} \nabla b_{k,i} \\ + \sum_{j=1}^{m} \left\{ \left[ \partial_t \mathcal{M}_{k,j} + \nabla \cdot (\mathcal{M}_{k,j} \mathbf{u}_k) - \frac{\partial \mathcal{L}_k}{\partial \xi_j} \right] \nabla \xi_j - (\nabla \nabla \xi_j) \mathbf{D}_{k,j} \right\} = 0. \quad (A.9) $$

We then use the fact that $(\nabla \nabla \xi_j) \mathbf{D}_{k,j} = \nabla \cdot (\nabla \xi_j \otimes \mathbf{D}_{k,j}) - (\nabla \cdot \mathbf{D}_{k,j}) \nabla \xi_j$ and rearrange the terms to obtain the generic momentum equations (2.5a) and (2.5b).

## A.2 Proof of Theorem 2.2

We provide a proof of the conservation law (2.8) for the Hamiltonian verified by the generic system. We derive the equations for the partial Hamiltonians (2.9). Equation (2.8) is then obtained by summation over the phases. Owing to the definition of the partial Hamiltonian $\mathcal{H}_k$ (2.7), we have

$$ m_k D_t^k \left( \frac{\mathcal{H}_k}{m_k} \right) = m_k D_t^k \left( \frac{\mathbf{K}_k^\top \mathbf{u}_k}{m_k} \right) + \sum_{j=1}^{m} m_k D_t^k \left( \frac{\mathcal{M}_{k,j} D_t^k \xi_j}{m_k} \right) - m_k D_t^k \left( \frac{\mathcal{L}_k}{m_k} \right). \quad (A.10) $$

We now derive a more suitable expression for $m_k D_t^k (\mathcal{L}_k/m_k) = D_t^k \mathcal{L}_k - \mathcal{L}_k D_t^k m_k/m_k$. Using the chain rule, we have

$$ D_t^k \mathcal{L}_k = m_k \frac{\partial \mathcal{L}_k}{\partial m_k} \frac{D_t^k m_k}{m_k} + \frac{\partial \mathcal{L}_k}{\partial \mathbf{u}_k} D_t^k \mathbf{u}_k + \sum_{i=1}^{n_{k'}} \frac{\partial \mathcal{L}_k}{\partial b_{k,i}} D_t^k b_{k,i} \\ + \sum_{j=1}^{m} \left\{ \frac{\partial \mathcal{L}_k}{\partial \xi_j} D_t^k \xi_j + \frac{\partial \mathcal{L}_k}{\partial (\nabla \xi_j)} D_t^k (\nabla \xi_j) + \frac{\partial \mathcal{L}_k}{\partial (D_t^k \xi_j)} D_t^k (D_t^k \xi_j) \right\} \quad (A.11) $$

Moreover, we have the following identities

$$ \frac{\partial \mathcal{L}_k}{\partial \mathbf{u}_k} D_t^k \mathbf{u}_k = m_k D_t^k \left( \frac{\mathbf{K}_k^\top \mathbf{u}_k}{m_k} \right) - m_k \mathbf{u}_k^\top D_t^k \left( \frac{\mathbf{K}_k}{m_k} \right), \quad (A.12) $$

$$ \frac{\partial \mathcal{L}_k}{\partial (D_t^k \xi_j)} D_t^k (D_t^k \xi_j) = m_k D_t^k \left( \frac{\mathcal{M}_{k,j} D_t^k \xi_j}{m_k} \right) - m_k (D_t^k \xi_j) D_t^k \left( \frac{\mathcal{M}_{k,j}}{m_k} \right). \quad (A.13) $$

Using these identities as well as the generic momentum equation for $\mathbf{K}_k$, it results that

$$ m_k D_t^k \left( \frac{\mathcal{L}_k}{m_k} \right) = m_k D_t^k \left( \frac{\mathbf{K}_k^\top \mathbf{u}_k}{m_k} \right) + \sum_{j=1}^{m} m_k D_t^k \left( \frac{\mathcal{M}_{k,j} D_t^k \xi_j}{m_k} \right) \underbrace{- \mathcal{L}_k^* (\nabla \cdot \mathbf{u}_k) - \mathbf{u}_k^\top \nabla \mathcal{L}_k^*}_{=T_1} \\ + \underbrace{\sum_{i=1}^{n_{k'}} \frac{\partial \mathcal{L}_k}{\partial b_{k',i}} D_t^k b_{k',i} + \sum_{i=1}^{n_k} \frac{\partial \mathcal{L}_{k'}}{\partial b_{k,i}} \mathbf{u}_k^\top \nabla b_{k,i} - \sum_{i=1}^{n_{k'}} \frac{\partial \mathcal{L}_k}{\partial b_{k',i}} \mathbf{u}_k^\top \nabla b_{k',i}}_{=T_2} \\ + \sum_{j=1}^{m} \left\{ \underbrace{\frac{\partial \mathcal{L}_k}{\partial \xi_j} D_t^k \xi_j - \frac{\partial \mathcal{L}_k}{\partial \xi_j} \mathbf{u}_k^\top \nabla \xi_j}_{=T_3} + m_k \mathbf{u}_k^\top \nabla \xi_j D_t^k \left( \frac{\mathcal{M}_{k,j}}{m_k} \right) - m_k D_t^k \xi_j D_t^k \left( \frac{\mathcal{M}_{k,j}}{m_k} \right)}_{=T_4} \\ \underbrace{\frac{\partial \mathcal{L}_k}{\partial (\nabla \xi_j)} D_t^k (\nabla \xi_j) - \mathbf{u}_k^\top \nabla \cdot (\nabla \xi_j \otimes \mathbf{D}_{k,j}) + \mathbf{u}_k^\top (\nabla \cdot \mathbf{D}_{k,j}) \nabla \xi_j}_{=T_5} \right\}. \quad (A.14) $$
Here, for convenience, we used to notation $\mathcal{L}_k^* = m_k \partial \mathcal{L}_k / \partial m_k - \mathcal{L}_k$. Terms $T_1$ to $T_5$ are then simplified into

$$ T_1 = -\nabla \cdot (\mathcal{L}_k^* \mathbf{u}_k), \qquad T_2 = \sum_{i=1}^{n_{k'}} \frac{\partial \mathcal{L}_k}{\partial b_{k',i}} \partial_t b_{k',i} - \sum_{i=1}^{n_k} \frac{\partial \mathcal{L}_{k'}}{\partial b_{k,i}} \partial_t b_{k,i}, \quad (A.15a) $$

$$ T_3 = \sum_{j=1}^{m} \frac{\partial \mathcal{L}_k}{\partial \xi_j} \partial_t \xi_j, \qquad T_4 = - \sum_{j=1}^{m} m_k D_t^k \left( \frac{\mathcal{M}_{k,j}}{m_k} \right) \partial_t \xi_j, \quad (A.15b) $$

$$ T_5 = \sum_{j=1}^{m} \{ \nabla \cdot (\mathbf{D}_{k,j} \partial_t \xi_j) - \nabla \cdot (\mathbf{D}_{k,j} \partial_t \xi_j) \}. \quad (A.15c) $$

Injecting these into (A.14), using the partial mass conservation and rearranging the various terms then yields

$$ m_k D_t^k \left( \frac{\mathcal{H}_k}{m_k} \right) - \nabla \cdot \left( \mathcal{L}_k^* \mathbf{u}_k + \sum_{j=1}^m \mathbf{D}_{k,j} \partial_t \xi_j \right) = \sum_{i=1}^{n_k} \frac{\partial \mathcal{L}_{k'}}{\partial b_{k,i}} \partial_t b_{k,i} - \sum_{i=1}^{n_{k'}} \frac{\partial \mathcal{L}_k}{\partial b_{k',i}} \partial_t b_{k',i} \\ + \sum_{j=1}^m \left\{ \partial_t \mathcal{M}_{k,j} + \nabla \cdot (\mathcal{M}_{k,j} \mathbf{u}_k + \mathbf{D}_{k,j}) - \frac{\partial \mathcal{L}_k}{\partial \xi_j} \right\} \partial_t \xi_j, \quad (A.16) $$

which, owing to the partial mass equation, is equivalent to (2.9).

## B Derivation of the new model in the Eulerian framework

We derive the new all-topology model (3.6) by Stationary Action Principle in the Eulerian framework. We consider the Lagrangian (3.9) and derive eqs. (3.3a) and (3.5a). This yields the isentropic version of the model, its energy form is then obtained through a straightforward change of variables to replace the phasic entropies $s_1$ and $s_2$ with the phasic total energy densities $m_1 E_1 = m_1 e_1 + m_1 u_1^2 / 2$ and $m_2 E_2 = m_2 e_2 + m_2 u_2^2 / 2$.

We consider the action defined through the Lagrangian (3.9), and assume independent and compactly supported perturbations for each variable. Using integrations by parts, one may show that for any pair of $\mathcal{C}^1$ functions $\phi, \psi$, where at least one has compact support in $\{t \in (t_0, t_1), \mathbf{x} \in \Omega_t\}$, we have

$$ \int_t \int_{\Omega_t} \phi D_t^k \psi \, d\mathbf{x} \, dt = - \int_t \int_{\Omega_t} \psi \left\{ \partial_t \phi + \nabla \cdot (\phi \mathbf{u}_k) \right\} \, d\mathbf{x} \, dt, \quad (B.1) $$

and an analogous result holds for material derivative and conservation at the mixture velocity $\mathbf{u}$. Using this result, the critical points of the Action satisfy the following equations, each obtained by one of the independent perturbations as specified,

$$ \delta\zeta: \quad D_t\alpha_1 = 0, \quad (B.2a) $$

$$ \delta\alpha_1: \quad \partial_t(\rho\zeta) + \nabla \cdot (\rho\zeta\mathbf{u}) = p_1 - p_2, \quad (B.2b) $$

$$ \delta\theta_k: \quad D_t^k s_k = 0, \quad (B.2c) $$

$$ \delta s_k: \quad \partial_t(m_k \theta_k) + \nabla \cdot (m_k \theta_k) = -m_k T_k, \quad (B.2d) $$

$$ \delta\lambda_{k,i}: \quad D_t^k \mu_{k,i} = 0, \quad (B.2e) $$

$$ \delta\mu_{k,i}: \quad \partial_t(m_k \lambda_{k,i}) + \nabla \cdot (m_k \lambda_{k,i} \mathbf{u}_k), \quad (B.2f) $$

$$ \delta\xi_k: \quad \partial_t m_k + \nabla \cdot (m_k \mathbf{u}_k) = 0, \quad (B.2g) $$

$$ \delta m_k: \quad D_t^k \xi_k - \theta_k D_t^k s_k - \zeta D_t \alpha_1 - \sum_i \lambda_{k,i} D_t^k \mu_{k,i} = \frac{1}{2} \mathbf{u}_k^2 - h_k, \quad (B.2h) $$

$$ \delta \mathbf{u}_k: \quad \mathbf{u}_k = \nabla \xi_k - \theta_k \nabla s_k - \zeta \nabla \alpha_1 - \sum_i \lambda_{k,i} \nabla \mu_{k,i}. \quad (B.2i) $$
In particular, eq. (B.2a) yields the volume fraction equation, while eq. (B.2g) yields the phasic mass conservation. Note that owing to the partial mass and total mass conservation laws, eqs. (B.2b), (B.2d) and (B.2f) may also be written

$$D_t \zeta = (p_1 - p_2)/\rho, \quad D_t^k \theta_k = -T_k, \quad D_t^k \lambda_{k,i} = 0, \qquad (B.3)$$

and so we obtain the phasic entropy equations, and by applying the gradient operator, we obtain the evolution equation for $v = \nabla \zeta$. At this stage, only the phasic momentum equations remain to be derived. Using the following identity,

$$D_t^k (\phi \nabla \psi) = \nabla (\phi D_t^k \psi) + (\nabla \psi) D_t^k \phi - (\nabla \phi) D_t^k \psi + (\nabla u_k) (\phi \nabla \psi), \qquad (B.4)$$

we obtain the relation

$$D_t^k \left( \nabla \xi_k - \theta_k \nabla s_k - \zeta \nabla \alpha_1 - \sum_i \lambda_{k,i} \nabla \mu_{k,i} \right) = \\ \nabla \left[ D_t^k \xi_k - \theta_k D_t^k s_k - \zeta D_t^k \alpha_1 - \sum_i \lambda_{k,i} D_t^k \mu_{k,i} \right] \\ + (\nabla \theta_k) D_t^k s_k - (\nabla s_k) D_t^k \theta_k + (\nabla \zeta) D_t^k \alpha_1 - (\nabla \alpha_1) D_t^k \zeta + \sum_i \left\{ (\nabla \lambda_{k,i}) D_t^k \mu_{k,i} - (\nabla \mu_{k,i}) D_t^k \lambda_{k,i} \right\} \\ - (\nabla u_k) \left( \nabla \xi_k - \theta_k \nabla s_k - \zeta \nabla \alpha_1 - \sum_i \lambda_{k,i} \nabla \mu_{k,i} \right). \qquad (B.5)$$

Using the decomposition (B.2i), and the evolutions equations given by (B.2h) and (B.3), the previous equation writes

$$D_t^k u_k + (\nabla u_k) u_k = \nabla \left( \frac{1}{2} u_k^2 - h_k \right) + T_k \nabla s_k + (\nabla \zeta) D_t^k \alpha_1 - (\nabla \alpha_1) D_t^k \zeta. \qquad (B.6)$$

We then use the thermodynamic identity $d h_k = T_k ds_k + \rho_k^{-1} dp_k$, and the following equations,

$$D_t^k \alpha_1 = D_t \alpha_1 + (\nabla \alpha_1)^\top (u_k - u) = (\nabla \alpha_1)^\top (u_k - u), \qquad (B.7)$$

$$D_t^k \zeta = D_t \zeta + (\nabla \zeta)^\top (u_k - u) = (p_1 - p_2)/\rho + (\nabla \zeta)^\top (u_k - u), \qquad (B.8)$$

to obtain the velocity equation

$$D_t^k u_k + \frac{1}{\rho_k} \nabla p_k + \frac{p_1 - p_2}{\rho} \nabla \alpha_1 = [v \otimes \nabla \alpha_1 - \nabla \alpha_1 \otimes v] (u_k - u), \qquad (B.9)$$

which, when multiplied by the partial mass $m_k$, yields the phasic momentum equation.

## References

[1] Loison A, Pichard T, Kokh S, Massot M. 2025 Small-Scale Interface Dynamic Modelling Based on the Geometric Method of Moments for a Two-Scale Two-Phase Flow Model with a Disperse Small Scale. *Journal of Fluid Mechanics* **1003**, A27.

[2] Frederix EMA, Dovizio D, Mathur A, Komen EMJ. 2021 All-Regime Two-Phase Flow Modeling Using a Novel Four-Field Large Interface Simulation Approach. *International Journal of Multiphase Flow* **145**, 103822.

[3] Loison A, Kokh S, Pichard T, Massot M. 2024 A Unified Two-Scale Gas-Liquid Multi-Fluid Model with Capillarity and Interface Regularization through a Mass Transfer between Scales. *International Journal of Multiphase Flow* **177**, 104857.
[4] Ishii M. 1975 *Thermo-Fluid Dynamic Theory of Two-phase Flow*. Paris: Eyrolles.

[5] Stuhmiller JH. 1977 The Influence of Interfacial Pressure Forces on the Character of Two-Phase Flow Model Equations. *International Journal of Multiphase Flow* **3**, 551–560.

[6] Ramshaw JD, Trapp JA. 1978 Characteristics, Stability, and Short-Wavelength Phenomena in Two-Phase Flow Equation Systems. *Nuclear Science and Engineering* **66**, 93–102.

[7] Vazquez-Gonzalez T, Llor A, Fochesato C. 2016 Ransom Test Results from Various Two-Fluid Schemes : Is Enforcing Hyperbolicity a Thermodynamically Consistent Option ?. *International Journal of Multiphase Flow* **81**, 104–112.

[8] Fox RO, Laurent F, Vié A. 2020 A Hyperbolic Two-Fluid Model for Compressible Flows with Arbitrary Material-Density Ratios. *Journal of Fluid Mechanics* **903**, A5.

[9] Ransom VH, Hicks DL. 1984 Hyperbolic Two-Pressure Models for Two-Phase Flow. *Journal of Computational Physics* **53**, 124–151.

[10] Baer MR, Nunziato JW. 1986 A Two-Phase Mixture Theory for the Deflagration-to-Detonation Transition (DDT) in Reactive Granular Materials. *International Journal of Multiphase Flow* **12**, 861–889.

[11] Chen Y, Glimm J, Sharp DH, Zhang Q. 1996 A Two-phase Flow Model of the Rayleigh–Taylor Mixing Zone. *Physics of Fluids* **8**, 816–825.

[12] Saurel R, Abgrall R. 1999 A Multiphase Godunov Method for Compressible Multifluid and Multiphase Flows. *Journal of Computational Physics* **150**, 425–467.

[13] Saurel R, Lemetayer O. 2001 A Multiphase Model for Compressible Flows with Interfaces, Shocks, Detonation Waves and Cavitation. *Journal of Fluid Mechanics* **431**, 239–271.

[14] Gallouët T, Hérard JM, Seguin N. 2004 Numerical Modeling of Two-Phase Flows Using the Two-Fluid Two-Pressure Approach. *Mathematical Models and Methods in Applied Sciences* **14**, 663–700.

[15] Tokareva SA, Toro EF. 2010 HLLC-type Riemann Solver for the Baer–Nunziato Equations of Compressible Two-Phase Flow. *Journal of Computational Physics* **229**, 3573–3604.

[16] Coquel F, Hérard JM, Saleh K, Seguin N. 2013 Two Properties of Two-Velocity Two-Pressure Models for Two-Phase Flows. *Communications in Mathematical Sciences* **12**, 593–600.

[17] Hantke M, Müller S. 2019 Closure Conditions for a One Temperature Non-Equilibrium Multi-Component Model of Baer-Nunziato Type. *ESAIM: Proceedings and Surveys* **66**, 42–60.

[18] Thanh MD, Vinh DX. 2022 The Riemann Problem for the Nonisentropic Baer–Nunziato Model of Two-Phase Flows. *Nonlinear Analysis: Real World Applications* **67**, 103623.

[19] Saurel R, Gavrilyuk S, Renaud F. 2003 A Multiphase Model with Internal Degrees of Freedom: Application to Shock-Bubble Interaction. *Journal of Fluid Mechanics* **495**, 283–321.

[20] Jin H, Glimm J, Sharp DH. 2006 Entropy of Averaging for Compressible Two-Pressure Two-Phase Flow Models. *Physics Letters A* **360**, 114–121.
[21] Perrier V, Gutiérrez E. 2021 Derivation and Closure of Baer and Nunziato Type Multiphase Models by Averaging a Simple Stochastic Model. *Multiscale Modeling & Simulation* **19**, 401–439.

[22] Guillemaud V. 2007 *Modélisation et simulation numérique des écoulements diphasiques par une approche bifluide à deux pressions*. PhD thesis Université de Provence - Aix-Marseille I.

[23] Lochon H, Daude F, Galon P, Hérard JM. 2016 Comparison of Two-Fluid Models on Steam-Water Transients. *ESAIM: Mathematical Modelling and Numerical Analysis* **50**, 1631–1657.

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

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

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

[27] Godunov S. 1961 An interesting class of quasilinear systems. In *Dokl. Akad. Nauk SSSR vol. 139 pp. 521–523*.

[28] Godunov S, Mikhailova T, Romenskii E. 1996 Systems of thermodynamically coordinated laws of conservation invariant under rotations. *Siberian Mathematical Journal* **37**, 690–705.

[29] Romensky E. 1998 Hyperbolic systems of thermodynamically compatible conservation laws in continuum mechanics. *Mathematical and computer modelling* **28**, 115–130.

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

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

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

[33] Haegeman W, Orlando G, Kokh S, Massot M. 2025 Modelling of Relative Velocity, Velocity Fluctuations and Their Interactions for Two-Fluid Models by Stationary Action Principle. *ESAIM: Proceedings and Surveys*. (submitted, preprint available on HAL hal-05224712v1).

[34] La Spina G, de' Michieli Vitturi M, Clarke AB. 2017 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.

[35] Gavrilyuk S. 2020 ‘Uncertainty’ Principle in Two Fluid–Mechanics. *ESAIM: Proceedings and Surveys* **69**, 47–55.

[36] Bretherton FP. 1970 A Note on Hamilton's Principle for Perfect Fluids. *Journal of Fluid Mechanics* **44**, 19–31.

[37] Salmon R. 1988 Hamiltonian Fluid Mechanics. *Annual Review of Fluid Mechanics* **20**, 225–256.
[38] Serre D. 1993 Sur le principe variationnel des équations de la mécanique des fluides parfaits. *ESAIM: Mathematical Modelling and Numerical Analysis* **27**, 739-758.

[39] Haegeman W, Kokh S, Massot M, Orlando G. 2025 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*. (In press, preprint available on HAL hal-05091807v2).

[40] Dhaouadi F, Gavrilyuk S. 2024 An Eulerian Hyperbolic Model for Heat Transfer Derived via Hamilton's Principle: Analytical and Numerical Study. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences* **480**, 20230440.

[41] Bedford A, Drumheller D. 1978 A Variational Theory of Immiscible Mixtures. *Archive for Rational Mechanics and Analysis* **68**, 37-51.

[42] Geurst JA. 1986 Variational Principles and Two-Fluid Hydrodynamics of Bubbly Liquid/Gas Mixtures. *Physica A: Statistical Mechanics and its Applications* **135**, 455-486.

[43] Gavrilyuk S, Saurel R. 2002 Mathematical and Numerical Modeling of Two-Phase Compressible Flows with Micro-Inertia. *Journal of Computational Physics* **175**, 326-360.

[44] Hérard JM. 2003 Numerical Modelling of Turbulent Two Phase Flows Using the Two Fluid Approach. In *16th AIAA Computational Fluid Dynamics Conference*, Fluid Dynamics and Co-located Conferences. American Institute of Aeronautics and Astronautics.

[45] Müller S, Hantke M, Richter P. 2016 Closure Conditions for Non-Equilibrium Multi-Component Models. *Continuum Mechanics and Thermodynamics* **28**, 1157-1189.

[46] Rubinow SI, Keller JB. 1961 The Transverse Force on a Spinning Sphere Moving in a Viscous Fluid. *Journal of Fluid Mechanics* **11**, 447-459.

[47] Crowe CT, Schwarzkopf JD, Sommerfeld M, Tsuji Y. 2011 *Multiphase Flows with Droplets and Particles*. Boca Raton: CRC Press 2 edition.

[48] Kapila AK, Menikoff R, Bdzil JB, Son SF, Stewart DS. 2001 Two-Phase Modeling of Deflagration-to-Detonation Transition in Granular Materials: Reduced Equations. *Physics of Fluids* **13**, 3002-3024.

[49] Petitpas F, Saurel R, Franquet E, Chinnayya A. 2009 Modelling Detonation Waves in Condensed Energetic Materials: Multiphase CJ Conditions and Multidimensional Computations. *Shock Waves* **19**, 377-401.

[50] Schmidmayer K, Cazé J, Petitpas F, Daniel É, Favrie N. 2023 Modelling Interactions between Waves and Diffused Interfaces. *International Journal for Numerical Methods in Fluids* **95**, 215-241.

[51] Hérard JM, Jomée G. 2023 Two Approaches to Compute Unsteady Compressible Two-Phase Flow Models with Stiff Relaxation Terms. *ESAIM: Mathematical Modelling and Numerical Analysis* **57**, 3537-3583.

[52] Labois M. 2008 *Modélisation des déséquilibres mécaniques pour les écoulements diphasiques : approches par relaxation et par modèle réduit*. PhD thesis Université de Provence - Aix-Marseille I. tel-00338818.

[53] Lund H. 2012 A Hierarchy of Relaxation Models for Two-Phase Flow. *SIAM Journal on Applied Mathematics* **72**, 1713-1741.
[54] Hantke M, Müller S, Grabowsky L. 2021 News on Baer–Nunziato-type Model at Pressure Equilibrium. *Continuum Mechanics and Thermodynamics* **33**, 1–22.

[55] Wood AB. 1930 *A Textbook of Sound*. London: G. Bell and Sons Ltd.

[56] Godlewski E, Raviart PA. 2021 *Numerical Approximation of Hyperbolic Systems of Conservation Laws* vol. 118 *Applied Mathematical Sciences*. New York, NY: Springer.

[57] Lhuillier D, Chang CH, Theofanous T. 2013 On the Quest for a Hyperbolic Effective-Field Model of Disperse Flows. *Journal of Fluid Mechanics* **731**.

[58] Menikoff R, Plohr BJ. 1989 The Riemann Problem for Fluid Flow of Real Materials. *Reviews of Modern Physics* **61**, 75–130.

[59] Thompson PA. 1971 A Fundamental Derivative in Gasdynamics. *The Physics of Fluids* **14**, 1843–1849.

[60] Friedrichs KO. 1954 Symmetric Hyperbolic Linear Differential Equations. *Communications on Pure and Applied Mathematics* **7**, 345–392.

[61] Kato T. 1975 The Cauchy Problem for Quasi-Linear Symmetric Hyperbolic Systems. *Archive for Rational Mechanics and Analysis* **58**, 181–205.

[62] Godunov SK. 2025 Interesting Class of Quasilinear Systems. *Journal of Computational Physics* **520**, 113521. Translation of the original manuscript published in Report of the USSR Academy of Sciences, 1961, Volume 139, No. 3, pp. 521–523.

[63] Friedrichs KO, Lax PD. 1971 Systems of Conservation Equations with a Convex Extension. *Proceedings of the National Academy of Sciences* **68**, 1686–1688.

[64] Forestier A, Gavrilyuk S. 2011 Criterion of Hyperbolicity for Non-Conservative Quasilinear Systems Admitting a Partially Convex Conservation Law. *Mathematical Methods in the Applied Sciences* **34**, 2148–2158.

[65] Dumbser M, Fambri F, Gaburro E, Reinarz A. 2020 On GLM Curl Cleaning for a First Order Reduction of the CCZ4 Formulation of the Einstein Field Equations. *Journal of Computational Physics* **404**, 109088.

[66] Godunov SK. 2025 Symmetric Form of the Magnetohydrodynamics Equations. *Journal of Computational Physics* **521**, 113523. Translation of the original manuscript published in Numerical methods of continuum mechanics, Academy of Sciences of the USSR, Siberian Branch, Computing Center, Volume 3, Issue 1, Novosibirsk, 1972.

[67] Powell KG, Roe PL, Linde TJ, Gombosi TI, De Zeeuw DL. 1999 A Solution-Adaptive Upwind Scheme for Ideal Magnetohydrodynamics. *Journal of Computational Physics* **154**, 284–309.