# Weak solution for granular model

Thierry Dubois, Laurent Chupin

► To cite this version:

Thierry Dubois, Laurent Chupin. Weak solution for granular model. 2025. hal-05080698v1

HAL Id: hal-05080698

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

Preprint submitted on 22 May 2025 (v1), last revised 3 Dec 2025 (v2)

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.
# Weak solution for granular model

Laurent Chupin * & Thierry Dubois *

## Abstract

This article is devoted to questions concerning the existence of solutions for partial differential equation problems modeling granular flows. The models studied take into account the complex threshold rheology of these flows, as well as the dilatance effects. It is the coupling of these two physical phenomena that ensures stability and the existence of dissipated energy. The key point of the article is to understand how this energy can ensure the existence of a weak solution. We first establish a complete result on a simplified model, then demonstrate how it can be extended to more general cases. This work represents a real breakthrough in the mathematical analysis of this type of models for complex flows.

**keywords:** granular model, weak existence, rheology, dilatance.

## 1 Introduction

**From grain to continuous model** Granular flow modeling is a major challenge in many contexts. Indeed, this type of flow is involved in a wide variety of fields, such as the earth sciences, with pyroclastic flows or dune movements; the food industry in the broadest sense of the term, from cereal silos to food design; medicine, and in particular the pharmaceutical industry, which makes extensive use of powders and capsules corresponding to granular media with diverse behaviors. The tricky aspect in this type of medium lies in its intermediate position between a fluid-like medium (such as gases, where the particles are too small to be individually described) and a solid medium, which requires understanding the interactions between a limited number of macroscopic particles.

With the advent of computational tools, we could be led to use these models even when the media comprise a very large number of particles (generally several billion). Currently this approach is hampered by our limits knowledge of the local interactions: at the grain level, the solid contact laws between multiple particles involve highly non-linear relations such as frictional or inelastic shocks, that remain difficult to fully understand. The well-known Discrete Element Method (see for instance [24]) attempts to take into account as much information as possible, but the lack of physical understanding means that not all the expected results can be obtained.

An alternative strategy is to define spatially averaged quantities to obtain continuous models, similar to what is done for fluid flows with the Navier-Stokes equations. One of the foundational works of this approach for granular flows is arguably that of the GdR MIDI, which highlighted the $\mu(I)$-rheology model as a basis for granular flow models (see, for instance, the collective work [21]).

**Issues and open questions** Since the $\mu(I)$-model is primarily derived from experimental considerations, numerous questions remain regarding its applicability in more general contexts those in which it was established (for instance, at larger scales). Furthermore, its implementation raises additional challenges to be used numerically. Indeed, to be able to numerically discretise a model and perform a computational implementation, it is highly desirable that it has a unique solution and that this solution is stable, in particular not too sensitive to perturbations in the data. Unfortunately, from this point of view, these continuous models often exhibit significant limitations including issues of instability (see the work of Barker and colleagues in [4, 5, 6]) as well as the absence of a comprehensive theory ensuring the existence of solutions. At first glance, these challenges appear to be closely linked to the non-linear and singular nature of the model. Despite these fundamental issues, the implementation of these models in a numerical framework is feasible,

<footnote>* Université Clermont Auvergne, CNRS, LMBP, F-63000 Clermont-Ferrand, France
Corresponding author: laurent.chupin@uca.fr</footnote>
but their reliability cannot be guaranteed. Indeed, instabilities have been frequently observed with classical formulations as reported in [15]. Several authors have presented results based on these models, but they have introduced assumptions that are not always physically relevant. This is the case, for example, in [1, 2, 16] where the authors examine the existence of a solution in a model with a pressure-dependent stress threshold. However, the pressure considered in these models differs from the total pressure which results in a simpler coupling. To address this limitation, one would need to incorporate more realistic models commonly used to describe granular flows (models capable of capturing phenomena such as normal stress differences, etc.). Such models have been developed and discussed in the review article [18] and it is this type of model that we aim to explore in the present work.

The major challenge is therefore to develop a model for dense granular flows that is both physically consistent and mathematically well-posed. A first step in this direction was taken in [14], where a stable and physically consistent model was proposed: the particularly remarkable aspect is that it is the coupling between rheology and dilatancy that ensures the model’s stability. However, the strong non-linearities of the model, along with the presence of thresholds in the stress expression raise important questions regarding the existence and uniqueness of solutions.

**Results and article outline** In this article, we propose a natural functional framework in which it is possible to establish the existence and uniqueness of a solution. Once again, it is the combination of several physical ingredients that guarantees these results. In the first part (Section 2), we describe the complete model and then focus on a simplified version that retains the two main physical features essential for ensuring the existence of a solution:

- the incorporation of threshold rheology, with a pressure dependant yield criterion, as in the $\mu(I)$-model;

- the consideration of dilatancy, as proposed in [27], in particular through a velocity field whose divergence, representing volume changes, is governed by a balanced between shear and pressure effects.

The second part form the core of the article and is devoted to the question of the existence of solutions for the proposed model. It is structured into several subsections. After introducing the mathematical notations in Subsection 3.1, the following subsection presents a key component of the main result: it provides the weak formulation of the problem studied. We observed in particular, that this formulation seems to eliminate certain threshold-related non-linearities. However, Proposition 1 shows that this formulation effectively captures the full physical model and allows us to recover these non-linear effects. The final subsections (Subsections 3.3, 3.4 and 3.5) are devoted to the proof of the existence result, with a particular focus on handling the remaining non-linear terms and establishing the positivity of pressure (a property that is physically intuitive, but rarely proven rigorously from a mathematical perspective). These results represent a real innovation in the mathematical analysis of this type of model for complex flows: prove that a granular model, including threshold rheology and dilatation, is mathematically well-posed thanks to these relevant physical laws.

In the final Section 4, we return to the full problem and show how the preceding analysis can be adapted accordingly a specific study of the $\mu(I)$-model highlights how it relates to the “simplified” version initially considered. This section concludes by pointing out that many questions remain open concerning this class of models.

## 2 Modeling dense granular flow

Guided by the following principles: make maximum use of physical principles (conservation laws) and phenomenological analyses (constitutive relations, equations of state); ensure the linear stability of the model, based on Barker’s stability conditions presented in [28]; and obtain an “energy” that dissipates over time, a complete dense granular flow model is proposed in [14]. This model couples four unknowns of the problem, namely the volume fraction scalar field $\phi$, the velocity vector field $\mathbf{u}$, the symmetric extra-stress tensor field $\boldsymbol{\sigma}$ and the pressure scalar field $p$. Depending on the rheology chosen for the granular part, the expression of the various terms differs, but in the “simplest” case (corresponding to a Druker-Prager rheology), the model
is written as follows

$$ \phi \rho_0 (\partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u}) + \nabla p - \operatorname{div}(2\nu(\phi, |\mathbf{Du}|)\mathbf{Du}) = \phi \rho_0 \mathbf{g} + \operatorname{div} \boldsymbol{\sigma}, \quad (1) $$

$$ \boldsymbol{\sigma} : \mathbf{Su} = 2\alpha(\phi)p|\mathbf{Su}|, \quad |\boldsymbol{\sigma}| \le \alpha(\phi)p, \quad \boldsymbol{\sigma} = \boldsymbol{\sigma}^{\mathrm{T}} \quad \text{and} \quad \operatorname{tr} \boldsymbol{\sigma} = 0, \quad (2) $$

$$ \operatorname{div} \mathbf{u} = 2\alpha(\phi)|\mathbf{Su}| - \beta(\phi)\sqrt{p}, \quad (3) $$

$$ \partial_t \phi + \operatorname{div}(\phi \mathbf{u}) = 0. \quad (4) $$

In this version g designates the gravity vector, $\rho_0$ is the grain density, $\phi_{\max}$ corresponds to the maximum volume fraction (of the order of 0.6 for experiment flows in laboratory, see [26]), $\alpha$ and $\beta$ are positive functions depending on $\phi$ and on several physical parameters such as average grain size, friction coefficient, etc. The function $\nu$ depends on the model choices: it reflects the possibly non-linear effect of viscous forces and depends on the strain-rate tensor corresponding to the symmetrical part of the velocity gradient: $\mathbf{Du} = \frac{1}{2}(\nabla \mathbf{u} + (\nabla \mathbf{u})^{\mathrm{T}})$. The rheology also makes appear the deviatoric strain-rate tensor $\mathbf{Su} = \mathbf{Du} - \frac{1}{3}(\operatorname{div} \mathbf{u})\mathbf{I}_3$. For all these tensors, the norm is defined by $|\mathbf{A}|^2 = \frac{1}{2}\mathbf{A} : \mathbf{A}$ where $\mathbf{A} : \mathbf{B} = \operatorname{tr}(\mathbf{A}^{\mathrm{T}} \mathbf{B})$.

The rheology introduced here corresponds to the Drucker-Prager rheology in the sense that the relations (2) describing the relationship between stress and deformation can be rewritten as follows

$$ \boldsymbol{\sigma} = \alpha(\phi)p \frac{\mathbf{Su}}{|\mathbf{Su}|} \quad \text{if } |\mathbf{Su}| \neq 0, \quad (5) $$

$$ |\boldsymbol{\sigma}| \le \alpha(\phi)p, \quad \boldsymbol{\sigma} = \boldsymbol{\sigma}^{\mathrm{T}} \quad \text{and} \quad \operatorname{tr} \boldsymbol{\sigma} = 0 \quad \text{if } |\mathbf{Su}| = 0. $$

With this formulation, the quantity $\alpha(\phi)$ can be seen as the threshold, often linked to the friction angle.

The other fundamental feature of the model (1)–(4) is the expression for the divergence of the velocity field, i.e. the equation (3). This is an interpretation of the law derived from the work of F. Radjaï and S. Roux in [27]. This law states that the local volume variation (i.e. $\operatorname{div} \mathbf{u}$) must be related to the volume fraction variation $\phi$. More precisely, $\operatorname{div} \mathbf{u}$ must be proportional to the deviation between $\phi$ and an equilibrium state $\phi_{\mathrm{eq}}$ depending on $|\mathbf{Su}|$ and on $p$, see Subsection 4.1 for more details.

**Remark 1.** As explained in [14], it is possible to choose more complex rheologies, such as the $\mu(I)$-rheology. In this case, however, the dilatation law (3) must be adapted. We will come back to the important case of $\mu(I)$-rheology in Section 4. However, for most granular models, the pressure $p$ must be positive (for instance, the very definition of the inertial number $I$ involves $\sqrt{p}$, see Section 4).

**Remark 2.** This model does not take into account the interstitial gas in which the grains move. In [14], this additional contribution is taken into account: essentially, it involves adding pressure forces $\nabla p_f$ due to the gas in (1), and providing an evolution law for this new pressure $p_f$. This evolution law is mainly a heat equation-type law, see [14] for more details. Note in particular that for this more complete model, as for the model presented here, energy is dissipated over time.

**Remark 3.** In [14], the viscosity $\nu$ is not taken into account, even though the authors mention this possibility in Remarks 3 and 4. In the case presented here, this viscosity is essential to ensure the regularity of the velocity field.

In the following section, in order to highlight the tricky and original points, we focus on a simpler case of this model. We will retain only the rheological and dilatational contributions, and we assume that the volume fraction is constant. The model to be studied is the following

$$ \phi_0 \rho_0 \partial_t \mathbf{u} + \nabla p - \operatorname{div}(2\nu_0 |\mathbf{Du}|\mathbf{Du}) = \mathbf{f} + \operatorname{div} \boldsymbol{\sigma}, \quad (6) $$

$$ \boldsymbol{\sigma} : \mathbf{Su} = 2\alpha_0 p |\mathbf{Su}|, \quad |\boldsymbol{\sigma}| \le \alpha_0 p, \quad \boldsymbol{\sigma} = \boldsymbol{\sigma}^{\mathrm{T}} \quad \text{and} \quad \operatorname{tr} \boldsymbol{\sigma} = 0, \quad (7) $$

$$ \operatorname{div} \mathbf{u} = 2\alpha_0 |\mathbf{Su}| - \beta_0 \sqrt{p}, \quad (8) $$

where $\mathbf{f}$ are the external forces, namely $\mathbf{f} = \phi_0 \rho_0 \mathbf{g}$. Note that this model is not a special case of the full model, in which $\phi$ would simply have been replaced by a constant $\phi_0$, since the velocity field $\mathbf{u}$ is not divergence free. We show in Subsection 4.2 that this “simplified” system (6)–(8) corresponds to the main order of the model (1)–(4) in a specific regime, in particular when $\phi$ is close to a constant itself close to $\phi_{\max}$. We propose in Section 4 an adaptation of these results to more complete cases, in particular taking into account the evolution of the volume fraction $\phi$. Finally, note the choice of non-linear viscosity. This choice is essentially technical, even if this type of viscosity may appear physically. In fact, as we shall see, it ensures sufficient velocity regularity to give meaning to the product $p|\mathbf{Su}|$.
# 3 Mathematical results

## 3.1 Mathematical framework and notations

To carry out the mathematical study, we assume that all constants are normalized. Of course, these assumptions do not qualitatively alter what follows. The model considered in this section is therefore written as

$$ \partial_t \mathbf{u} + \nabla p - \mathrm{div}(2|\mathbf{Du}|\mathbf{Du}) = \mathbf{f} + \mathrm{div} \, \boldsymbol{\sigma}, \quad (9) $$

$$ \boldsymbol{\sigma} : \mathbf{Su} = 2p|\mathbf{Su}|, \quad |\boldsymbol{\sigma}| \le p, \quad \boldsymbol{\sigma} = \boldsymbol{\sigma}^T \quad \text{and} \quad \mathrm{tr} \, \boldsymbol{\sigma} = 0, \quad (10) $$

$$ \mathrm{div} \, \mathbf{u} = 2|\mathbf{Su}| - \sqrt{p}. \quad (11) $$

Given a time $T > 0$ and a bounded open domain $\Omega \subset \mathbb{R}^3$ whose boundary will be assumed to be of class $C^2$, we consider the problem (9)-(11) on $(0, T) \times \Omega$. It is completed by the following initial and boundary conditions:

$$ \mathbf{u}|_{t=0} = \mathbf{u}_{\mathrm{init}} \quad \text{and} \quad \mathbf{u}|_{\partial \Omega} = \mathbf{0}. \quad (12) $$

**Notations** - In order to define the notion of weak solution and to formulate the results, we need to fix the notations. The symbol $\| \cdot \|_q$ for $1 \le q \le +\infty$ stands for the $L^q$-norm in the usual Lebesgue space $L^q(\Omega)$ while $\| \cdot \|_{s,q}$ refers to the norm of the Sobolev space $W^{s,q}(\Omega)$, $s \in \mathbb{R}$ and $1 \le q \le +\infty$.

We will make repeated use of the Sobolev space $W_0^{1,3}(\Omega)$ of functions of $W^{1,3}(\Omega)$ which are zero at the boundary, as well as its dual $W^{-1,3/2}(\Omega)$. The associated duality bracket will be denoted $\langle \cdot, \cdot \rangle$.

Finally, for a Banach space $X$, we denote the relevant Bochner space by $L^q(0, T; X)$. The associated norm will be explicitly denoted $\| \cdot \|_{L^q(0,T;X)}$.

**Remark 4.** — As explained in the previous section, the relation (11) is derived from physical observations reflecting a competition between shear and pressure effects on volume changes. Nevertheless, we can also see a fundamental mathematical ingredient that helps us understand the link between this relation (11) and the rheology given by (10). Indeed, like the role of pressure in the Stokes equations describing Newtonian fluid flow, the pressure $p$ introduced in the equations (9)-(10) seems to be linked to a Lagrange multiplier associated with the constraint relation (11).

As an example, consider the minimization problem $\mathbf{u} = \underset{\mathbf{v}}{\operatorname{argmin}} \{f(\mathbf{v}); g(\mathbf{v}) = \mathbf{0}\}$ where

$$ f : \mathbf{v} \in W_0^{1,2}(\Omega) \mapsto \frac{1}{2} \int_{\Omega} |\nabla \mathbf{v}|^2 \in \mathbb{R} $$

$$ \text{and} \quad g : \mathbf{v} \in W_0^{1,2}(\Omega) \mapsto \mathrm{div} \, \mathbf{v} - 2|\mathbf{Sv}| - q \in L^2(\Omega), $$

where the source term $q$ is given in $L^2(\Omega)$. Note that the applications $f$ and $g$ are differentiable in any non-canceling function $\mathbf{v} \in W_0^{1,2}(\Omega)$ (if $\mathbf{v}$ cancels, we should probably consider the sub-differential of $g$, as proposed by Beck [7]):

$$ df(\mathbf{v}) : \varphi \in W_0^{1,2}(\Omega) \mapsto \int_{\Omega} \nabla \mathbf{v} : \nabla \varphi \in \mathbb{R} $$

$$ \text{and} \quad dg(\mathbf{v}) : \varphi \in W_0^{1,2}(\Omega) \mapsto \mathrm{div} \, \varphi - \frac{\mathbf{Sv}}{|\mathbf{Sv}|} : \mathbf{S}\varphi \in L^2(\Omega). $$

If a solution $\mathbf{u}$, minimizing $f$ under the constraint $g(\mathbf{u}) = \mathbf{0}$, exists and does not cancel, then according to the Lagrange multiplier Theorem, there exists a linear application $\tilde{p} : L^2(\Omega) \to \mathbb{R}$ such that $df(\mathbf{u}) = \tilde{p} \circ dg(\mathbf{u})$. According to the Riesz representation Theorem, the linear application $\tilde{p}$ can be represented by the scalar product in $L^2(\Omega)$: there exists $p \in L^2(\Omega)$ such that for all $\psi \in L^2(\Omega)$ we have $\tilde{p}(\psi) = \int_{\Omega} p\psi$. Thus, using the expressions for the differentials of the applications $f$ and $g$, we deduce that

$$ \forall \varphi \in W_0^{1,2}(\Omega) \quad \int_{\Omega} \nabla \mathbf{u} : \nabla \varphi = \int_{\Omega} p(\mathrm{div} \, \varphi - \frac{\mathbf{Sv}}{|\mathbf{Sv}|} : \mathbf{S}\varphi) $$

which indicates that, in the sense of distributions, we have $-\Delta \mathbf{u} + \nabla p = \mathrm{div} \, \boldsymbol{\sigma}$ with $\boldsymbol{\sigma} : \mathbf{Su} = 2p|\mathbf{Su}|$.

The approach described in Remark 4, which allows pressure to be seen as a Lagrange multiplier, will not be used in this article, but could perhaps provide some interesting results. We propose here an original formulation in terms of velocity-pressure, in which we emphasize the energy aspect.
## 3.2 Weak formulation and statement of result

An important point is the definition of the notion of solution of the system of equations (9)-(12).

**Definition 1** (weak solution). Let $\mathbf{u}_{\text{init}} \in L^2(\Omega)$ and $\mathbf{f} \in L^{\frac{3}{2}}(0, T; W^{-1, \frac{3}{2}}(\Omega))$.

We say that $(\mathbf{u}, p, \boldsymbol{\sigma})$ is a weak solution of (9)-(12) if

$$ \mathbf{u} \in L^{\infty}(0, T; L^2(\Omega)) \cap L^3(0, T; W_0^{1,3}(\Omega)), \quad \partial_t \mathbf{u} \in L^{\frac{3}{2}}(0, T; W^{-1, \frac{3}{2}}(\Omega)), $$

$$ p \in L^{\frac{3}{2}}((0, T) \times \Omega) \quad \text{with } p \ge 0 \quad \text{a.e.} $$

$$ \boldsymbol{\sigma} \in L^{\frac{3}{2}}((0, T) \times \Omega) \quad \text{with } |\boldsymbol{\sigma}| \le p, \quad \boldsymbol{\sigma} = \boldsymbol{\sigma}^{\top} \quad \text{and} \quad \text{tr} \boldsymbol{\sigma} = 0 \quad \text{a.e.} $$

and for all $\varphi \in L^3(0, T; W_0^{1,3}(\Omega))$, for all $\psi \in L^{\frac{3}{2}}((0, T) \times \Omega)$ we have

$$ \langle \partial_t \mathbf{u}, \varphi \rangle - \int_{\Omega} p \operatorname{div} \varphi + \int_{\Omega} 2|\mathbf{D}\mathbf{u}|\mathbf{D}\mathbf{u} : \mathbf{D}\varphi + \int_{\Omega} \boldsymbol{\sigma} : \mathbf{S}\varphi = \langle \mathbf{f}, \varphi \rangle, \quad (13) $$

$$ \int_{\Omega} \psi \operatorname{div} \mathbf{u} - \int_{\Omega} 2\psi|\mathbf{S}\mathbf{u}| + \int_{\Omega} \psi\sqrt{p} = 0, \quad (14) $$

$$ \frac{1}{2} \int_{\Omega} |\mathbf{u}|^2 + 4 \int_{0}^{T} \int_{\Omega} |\mathbf{D}\mathbf{u}|^3 + \int_{0}^{T} \int_{\Omega} |p|^{\frac{3}{2}} \le \frac{1}{2} \int_{\Omega} |\mathbf{u}_{\text{init}}|^2 + \int_{0}^{T} \langle \mathbf{f}, \mathbf{u} \rangle, \quad (15) $$

and the initial condition $\mathbf{u}|_{t=0} = \mathbf{u}_{\text{init}}$ holds in $L^2(\Omega)$.

Relationships (13) and (14) naturally translate equations (9) and (11) into their weak form. However it is not completely clear that threshold rheology as expressed in relations (10) is fully taken into account in this notion of solution.

**Proposition 1.** If $(\mathbf{u}, p, \boldsymbol{\sigma})$ is a weak solution of (9)-(12) then the relations (10) are satisfied almost everywhere.

**Proof** - Let $(\mathbf{u}, p, \boldsymbol{\sigma})$ be a weak solution of (9)-(12). Taking $\varphi = \mathbf{u}$ as test function in (13), $\psi = p$ in (14) and adding the results, we deduce

$$ \frac{1}{2} \frac{d}{dt} \int_{\Omega} |\mathbf{u}|^2 + \int_{\Omega} 4|\mathbf{D}\mathbf{u}|^3 + \int_{\Omega} |p|^{\frac{3}{2}} + \int_{\Omega} \left( \boldsymbol{\sigma} : \mathbf{S}\mathbf{u} - 2p|\mathbf{S}\mathbf{u}| \right) = \langle \mathbf{f}, \mathbf{u} \rangle. $$

Integrating with respect to time and comparing with (15), we get

$$ \int_{0}^{T} \int_{\Omega} \left( \boldsymbol{\sigma} : \mathbf{S}\mathbf{u} - 2p|\mathbf{S}\mathbf{u}| \right) \ge 0. \quad (16) $$

But, from the Cauchy-Schwarz inequality $\boldsymbol{\sigma} : \mathbf{S}\mathbf{u} \le 2|\boldsymbol{\sigma}||\mathbf{S}\mathbf{u}|$ and using the fact that $|\boldsymbol{\sigma}| \le p$, we have $\boldsymbol{\sigma} : \mathbf{S}\mathbf{u} - 2p|\mathbf{S}\mathbf{u}| \le 0$. Regarding (16) that implies $\boldsymbol{\sigma} : \mathbf{S}\mathbf{u} = 2p|\mathbf{S}\mathbf{u}|$ almost everywhere. ■

The main result is

**Theorem 1.** If $\mathbf{u}_{\text{init}} \in L^2(\Omega)$ and $\mathbf{f} \in L^{\frac{3}{2}}(0, T; W^{-1, \frac{3}{2}}(\Omega))$ then there exists a weak solution $(\mathbf{u}, p, \boldsymbol{\sigma})$ to the problem (9)-(12). The velocity $\mathbf{u}$ and the pressure $p$ are uniquely determined.

## 3.3 Approximations system

In order to obtain a weak solution, we start by constructing a solution to an approximate problem. To do this, consider two sequences $(\mathbf{f}_{\varepsilon})_{\varepsilon>0}$ and $(\mathbf{u}_{\text{init},\varepsilon})_{\varepsilon>0}$ of regular functions respectively converging to $\mathbf{f}$ in $L^{\frac{3}{2}}(0, T; W^{-1, \frac{3}{2}}(\Omega))$ and $\mathbf{u}_{\text{init}}$ in $L^2(\Omega)$. For $\varepsilon > 0$ we consider the following problem

$$ \partial_t \mathbf{u}_\varepsilon + \nabla p_\varepsilon - \operatorname{div}(2|\mathbf{D}\mathbf{u}_\varepsilon|\mathbf{D}\mathbf{u}_\varepsilon) = \mathbf{f}_\varepsilon + \operatorname{div} \boldsymbol{\sigma}_\varepsilon, \quad (17) $$

$$ \boldsymbol{\sigma}_\varepsilon = p_\varepsilon \frac{\mathbf{S}\mathbf{u}_\varepsilon}{|\mathbf{S}\mathbf{u}_\varepsilon| + \varepsilon}, \quad (18) $$

$$ \operatorname{div} \mathbf{u}_\varepsilon = 2 \frac{|\mathbf{S}\mathbf{u}_\varepsilon|^2}{|\mathbf{S}\mathbf{u}_\varepsilon| + \varepsilon} - V_\varepsilon(p_\varepsilon) - \varepsilon(\partial_t p_\varepsilon - \Delta p_\varepsilon). \quad (19) $$
with the following initial and boundary conditions:

$$ \mathbf{u}_{\varepsilon}|_{t=0} = \mathbf{u}_{\text{init},\varepsilon}, \quad \mathbf{u}_{\varepsilon}|_{\partial\Omega} = \mathbf{0}, \quad p_{\varepsilon}|_{t=0} = 0 \quad \text{and} \quad p_{\varepsilon}|_{\partial\Omega} = 0. \qquad (20) $$

The function $V_{\varepsilon}$ is an approximation of the square-root function, which allows us to define pressure $p_{\varepsilon}$ for all $\varepsilon > 0$, and to show that we obtain a positive pressure at the limit $\varepsilon \to 0$. More precisely, $V_{\varepsilon}$ is a concave function of class $C^1$, defined on $\mathbb{R}$ by

$$ V_{\varepsilon}(x) = \begin{cases} \sqrt{x + (\varepsilon/2)^2} - \varepsilon/2 & \text{if } x > 0, \\ x/\varepsilon & \text{if } x \le 0. \end{cases} \qquad (21) $$

The following two graphs illustrate how the square root function and the Heaviside function have been approximated by smooth functions.

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

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

**Proposition 2.** For $\varepsilon > 0$, there exists a unique weak solution $(\mathbf{u}_{\varepsilon}, p_{\varepsilon}, \boldsymbol{\sigma}_{\varepsilon})$ to equations (17)–(20), that is such that

$$ \begin{align*} \mathbf{u}_{\varepsilon} &\in L^{\infty}(0, T; L^2(\Omega)) \cap L^3(0, T; W_0^{1,3}(\Omega)), & \partial_t \mathbf{u}_{\varepsilon} &\in L^{\frac{3}{2}}(0, T; W^{-1, \frac{3}{2}}(\Omega)), \\ p_{\varepsilon} &\in L^{\infty}(0, T; L^2(\Omega)) \cap L^2(0, T; W_0^{1,2}(\Omega)), & \partial_t p_{\varepsilon} &\in L^2(0, T; W^{-1,2}(\Omega)), \\ \boldsymbol{\sigma}_{\varepsilon} &\in L^{\frac{3}{2}}((0, T) \times \Omega), & & & \end{align*} $$

with $\boldsymbol{\sigma}_{\varepsilon} = p_{\varepsilon} \frac{\mathbf{S}u_{\varepsilon}}{|\mathbf{S}u_{\varepsilon}| + \varepsilon}$ and for all $\varphi \in L^3(0, T; W_0^{1,3}(\Omega))$, for all $\psi \in L^2(0, T; W_0^{1,2}(\Omega))$

$$ \langle \partial_t \mathbf{u}_{\varepsilon}, \varphi \rangle - \int_{\Omega} p_{\varepsilon} \operatorname{div} \varphi + \int_{\Omega} 2|\mathbf{D}\mathbf{u}_{\varepsilon}|\mathbf{D}\mathbf{u}_{\varepsilon} : \mathbf{D}\varphi + \int_{\Omega} \boldsymbol{\sigma}_{\varepsilon} : \mathbf{S}\varphi = \langle \mathbf{f}_{\varepsilon}, \varphi \rangle, \qquad (22) $$

$$ \langle \varepsilon \partial_t p_{\varepsilon}, \psi \rangle + \int_{\Omega} \varepsilon \nabla p_{\varepsilon} \cdot \nabla \psi + \int_{\Omega} \psi \operatorname{div} \mathbf{u}_{\varepsilon} - \int_{\Omega} 2\psi \frac{|\mathbf{S}u_{\varepsilon}|^2}{|\mathbf{S}u_{\varepsilon}| + \varepsilon} + \int_{\Omega} \psi V_{\varepsilon}(p_{\varepsilon}) = 0. \qquad (23) $$

The initial conditions are satisfied in $L^2(\Omega): \mathbf{u}_{\varepsilon}|_{t=0} = \mathbf{u}_{\text{init},\varepsilon}$ and $p_{\varepsilon}|_{t=0} = 0$.

**Proof** - Adding $\partial_t p_{\varepsilon} - \Delta p_{\varepsilon}$ into the equation (19) allows us to define, for each $\varepsilon > 0$, the pressure $p_{\varepsilon}$ as the solution to a heat-type equation. The explicit expression of $\boldsymbol{\sigma}_{\varepsilon}$ given by (18) is a usual regularized version of the condition (10). This type of regularization is very common in this kind of problem, see for example [2]. The existence of a unique solution to the complete system (17)–(19) is then classical (for example, using a Galerkin method, see [2]). In practice, the Galerkin method applies as soon as we have an energy estimate for the solution. We therefore show below how to derive such an estimate.

We choose $\varphi = \mathbf{u}_{\varepsilon}$ as test function in (22), $\psi = p_{\varepsilon}$ in (23), we integrate with respect to time and add the results. Note, in particular that, based on the expression for $\boldsymbol{\sigma}_{\varepsilon}$, we have

$$ \boldsymbol{\sigma}_{\varepsilon} : \mathbf{S}u_{\varepsilon} = 2p_{\varepsilon} \frac{|\mathbf{S}u_{\varepsilon}|^2}{|\mathbf{S}u_{\varepsilon}| + \varepsilon}. $$

We also note that, for $\mathbf{u}_{\varepsilon} \in L^3(0, T; W_0^{1,3}(\Omega))$ and $\partial_t \mathbf{u}_{\varepsilon} \in L^{\frac{3}{2}}(0, T; W^{-1, \frac{3}{2}}(\Omega))$, we know that (see [11, p.99])

$$ \langle \partial_t \mathbf{u}_{\varepsilon}, \mathbf{u}_{\varepsilon} \rangle = \frac{1}{2} \frac{d}{dt} \|\mathbf{u}_{\varepsilon}\|_{2}^{2}. $$
The same reasoning applied to the time derivative of the pressure term finally implies

$$ \frac{1}{2} \frac{d}{dt} \|\mathbf{u}_{\varepsilon}\|_{2}^{2} + \frac{\varepsilon}{2} \frac{d}{dt} \|p_{\varepsilon}\|_{2}^{2} + 4 \|\mathbf{D} \mathbf{u}_{\varepsilon}\|_{3}^{3} + \varepsilon \|\nabla p_{\varepsilon}\|_{2}^{2} + \int_{\Omega} p_{\varepsilon} V_{\varepsilon}(p_{\varepsilon}) = \langle \mathbf{f}_{\varepsilon}, \mathbf{u}_{\varepsilon} \rangle. \quad (24) $$

This estimate allows us to demonstrate the existence of a solution to equations (17)–(20) in the classical way. Note in particular that the choice of the function $V_{\varepsilon}$ ensures that for all $x \in \mathbb{R}$, we have $xV_{\varepsilon}(x) \ge 0$. The proof of uniqueness is also standard, and we refer the reader to Subsection 3.7 for a similar reasoning leading to the result. ■

## 3.4 Uniform estimates

To take the limit as $\varepsilon \to 0$ in the system (17)–(20), we will establish uniform bounds with respect to $\varepsilon$ for the unknowns $\mathbf{u}_{\varepsilon}$, $p_{\varepsilon}$ and $\boldsymbol{\sigma}_{\varepsilon}$:

**Proposition 3.** The following quantities are bounded independently of $\varepsilon > 0$:

$$ \begin{aligned} & \|\mathbf{u}_{\varepsilon}\|_{L^{\infty}(0,T;L^2(\Omega))}, \quad \|\nabla \mathbf{u}_{\varepsilon}\|_{L^3((0,T)\times\Omega)}, \quad \|\partial_t \mathbf{u}_{\varepsilon}\|_{L^{\frac{3}{2}}(0,T;W^{-1,\frac{3}{2}}(\Omega))}, \quad \|\boldsymbol{\sigma}_{\varepsilon}\|_{L^{\frac{3}{2}}((0,T)\times\Omega)}, \\ & \|p_{\varepsilon}\|_{L^{\frac{3}{2}}((0,T)\times\Omega)}, \quad \|\varepsilon^{\frac{1}{2}} p_{\varepsilon}\|_{L^{\infty}(0,T;L^2(\Omega))}, \quad \|\varepsilon^{\frac{1}{2}} \nabla p_{\varepsilon}\|_{L^2((0,T)\times\Omega)}, \quad \|\varepsilon \partial_t p_{\varepsilon}\|_{L^2(0,T;W^{-1,2}(\Omega))}. \end{aligned} $$

**Proof - Step 1: using energy estimate**

We first estimate the term in the right-hand side of equation (24). Thus, using successively the duality between $W^{-1,\frac{3}{2}}(\Omega)$ and $W_0^{1,3}(\Omega)$, the Korn's inequality $\|\mathbf{u}\|_{1,3} \le K_3 \|\mathbf{D}\mathbf{u}\|_3$, see [20, Chapter V, Theorem 1.10], and the Young's inequality $ab \le \frac{2}{9}a^{\frac{3}{2}} + 3b^3$, we obtain

$$ \begin{aligned} \langle \mathbf{f}_{\varepsilon}, \mathbf{u}_{\varepsilon} \rangle &\le \|\mathbf{f}_{\varepsilon}\|_{-1, \frac{3}{2}} \|\mathbf{u}_{\varepsilon}\|_{1,3} \\ &\le K_3 \|\mathbf{f}_{\varepsilon}\|_{-1, \frac{3}{2}} \|\mathbf{D}\mathbf{u}_{\varepsilon}\|_3 \\ &\le \frac{2}{9} K_3^{\frac{3}{2}} \|\mathbf{f}_{\varepsilon}\|_{-1, \frac{3}{2}}^{\frac{3}{2}} + 3 \|\mathbf{D}\mathbf{u}_{\varepsilon}\|_3^3. \end{aligned} $$

Finally, the equality of energy (24) implies the following estimate:

$$ \frac{d}{dt} \|\mathbf{u}_{\varepsilon}\|_{2}^{2} + \varepsilon \frac{d}{dt} \|p_{\varepsilon}\|_{2}^{2} + 2 \|\mathbf{D}\mathbf{u}_{\varepsilon}\|_{3}^{3} + 2\varepsilon \|\nabla p_{\varepsilon}\|_{2}^{2} + 2 \int_{\Omega} p_{\varepsilon} V_{\varepsilon}(p_{\varepsilon}) \le \frac{4}{9} K_3^{\frac{3}{2}} \|\mathbf{f}_{\varepsilon}\|_{-1, \frac{3}{2}}^{\frac{3}{2}}. $$

Integrating in time and using the convergence of $\mathbf{f}_{\varepsilon}$ to $\mathbf{f}$ in $L^{\frac{3}{2}}(0, T; W^{-1, \frac{3}{2}}(\Omega))$, we deduce

$$ \begin{aligned} & \sup_{(0,T)} \|\mathbf{u}_{\varepsilon}\|_{2}^{2} + \varepsilon \sup_{(0,T)} \|p_{\varepsilon}\|_{2}^{2} + 2 \int_{0}^{T} \|\mathbf{D}\mathbf{u}_{\varepsilon}\|_{3}^{3} + 2\varepsilon \int_{0}^{T} \|\nabla p_{\varepsilon}\|_{2}^{2} + 2 \int_{0}^{T} \int_{\Omega} p_{\varepsilon} V_{\varepsilon}(p_{\varepsilon}) \\ & \le \|\mathbf{u}_{\text{init}}\|_{2}^{2} + K_3^{\frac{3}{2}} \int_{0}^{T} \|\mathbf{f}\|_{-1, \frac{3}{2}}^{\frac{3}{2}}. \end{aligned} \quad (25) $$

Noting that, by definition of the function $V_{\varepsilon}$, see (21), the quantity $p_{\varepsilon}V_{\varepsilon}(p_{\varepsilon})$ is always non-negative we therefore have estimates of the following quantities, independently of $\varepsilon$:

$$ \|\mathbf{u}_{\varepsilon}\|_{L^{\infty}(0,T;L^2(\Omega))}, \quad \|\varepsilon^{\frac{1}{2}} p_{\varepsilon}\|_{L^{\infty}(0,T;L^2(\Omega))}, \quad \|\mathbf{D}\mathbf{u}_{\varepsilon}\|_{L^3((0,T)\times\Omega)} \quad \text{and} \quad \|\varepsilon^{\frac{1}{2}} \nabla p_{\varepsilon}\|_{L^2((0,T)\times\Omega)}. $$

In addition, the estimate of the quantity $p_{\varepsilon}V_{\varepsilon}(p_{\varepsilon})$ provides a bound on the pressure $p_{\varepsilon}$. Indeed, we have

$$ \int_{0}^{T} \int_{\Omega} |p_{\varepsilon}|^{\frac{3}{2}} = \underbrace{\iint_{0 \le p_{\varepsilon} \le 1}}_{A} |p_{\varepsilon}|^{\frac{3}{2}} + \underbrace{\iint_{p_{\varepsilon} < 0}}_{B} |p_{\varepsilon}|^{\frac{3}{2}} + \underbrace{\iint_{p_{\varepsilon} > 1}}_{C} |p_{\varepsilon}|^{\frac{3}{2}}. $$

We clearly have $A \le T|\Omega|$. Moreover, due to Hölder's inequality, we obtain

$$ B \le \left( \iint_{p_{\varepsilon} < 0} \varepsilon^3 \right)^{\frac{1}{4}} \left( \iint_{p_{\varepsilon} < 0} \frac{1}{\varepsilon} |p_{\varepsilon}|^2 \right)^{\frac{3}{4}} \le (T|\Omega|\varepsilon^3)^{\frac{1}{4}} \left( \int_{0}^{T} \int_{\Omega} p_{\varepsilon} V_{\varepsilon}(p_{\varepsilon}) \right)^{\frac{3}{4}}. $$
Using the fact that, for $x \ge 1$ and $\varepsilon \le 2$ we have $x^{3/2} \le \frac{1}{\sqrt{2-1}}xV_\varepsilon(x)$, this implies that

$$C \le \frac{1}{\sqrt{2}-1} \int_0^T \int_\Omega p_\varepsilon V_\varepsilon(p_\varepsilon).$$

These three estimates on A, B and C combined with the energy estimate (25) provide a bound on $\|p_\varepsilon\|_{L^{3/2}((0,T) \times \Omega)}$ that is independent of $\varepsilon$, as long as $\varepsilon \le 2$).

The estimate on the stress $\sigma_\varepsilon$ follows directly from that of the pressure, since expression (18) yields the inequality $|\sigma_\varepsilon| \le |p_\varepsilon|$.

### Step 2: Bounds for times derivatives

Let us go back to the weak formulation (22). More precisely, using Hölder's inequality and the duality $W^{-1, 3/2}(\Omega) \leftrightarrow W_0^{1,3}(\Omega)$, we know that for all $\varphi \in L^3(0, T; W_0^{1,3}(\Omega))$ we have

$$\langle \partial_t \mathbf{u}_\varepsilon, \varphi \rangle \le \|p_\varepsilon\|_{3/2} \|\mathrm{div} \, \varphi\|_3 + 2 \|\mathbf{D} \mathbf{u}_\varepsilon\|_{3/2}^2 \|\mathbf{D} \varphi\|_3 + \|\sigma_\varepsilon\|_{3/2} \|\mathbf{S} \varphi\|_3 + \|\mathbf{f}_\varepsilon\|_{-1, 3/2} \|\varphi\|_{1,3}.$$

Integrating in time and applying Hölder's inequalities (with respect to integration in time), we derive the following estimate

$$\begin{aligned} \left| \int_0^T \langle \partial_t \mathbf{u}_\varepsilon, \varphi \rangle \right| \le & \left( \int_0^T \|p_\varepsilon\|_{3/2}^{3/2} \right)^{3/2} \left( \int_0^T \|\mathrm{div} \, \varphi\|_3^3 \right)^{1/3} + 2 \left( \int_0^T \|\mathbf{D} \mathbf{u}_\varepsilon\|_3^3 \right)^{3/2} \left( \int_0^T \|\mathbf{D} \varphi\|_3^3 \right)^{1/3} \\ & + \left( \int_0^T \|\sigma_\varepsilon\|_{3/2}^{3/2} \right)^{3/2} \left( \int_0^T \|\mathbf{S} \varphi\|_3^3 \right)^{1/3} + \left( \int_0^T \|\mathbf{f}_\varepsilon\|_{-1, 3/2}^{3/2} \right)^{3/2} \left( \int_0^T \|\varphi\|_{1,3}^3 \right)^{1/3}. \end{aligned}$$

By virtue of the estimates (25), we deduce that there exists a constant $C_1$, independent of $\varepsilon$, such that

$$\left| \int_0^T \langle \partial_t \mathbf{u}_\varepsilon, \varphi \rangle \right| \le C_1 \left( \int_0^T \|\varphi\|_{1,3}^3 \right)^{1/3} = C_1 \|\varphi\|_{L^3(0,T;W_0^{1,3}(\Omega))}.$$

Thus, by the definition of the dual norm, we obtain $\|\partial_t \mathbf{u}_\varepsilon\|_{L^{3/2}(0,T;W^{-1,3/2}(\Omega))} \le C_1$.

Similarly, we use (19) to estimate $\partial_t p_\varepsilon$. For all $\psi \in L^2(0, T; W_0^{1,2}(\Omega))$, we have

$$\langle \varepsilon \partial_t p_\varepsilon, \psi \rangle \le \varepsilon \|\nabla p_\varepsilon\|_2 \|\nabla \psi\|_2 + \|\psi\|_{3/2} \|\mathrm{div} \, \mathbf{u}_\varepsilon\|_3 + 2 \|\psi\|_{3/2} \|\mathbf{S} \mathbf{u}_\varepsilon\|_3 + \|\psi\|_{3/2} \|V_\varepsilon(p_\varepsilon)\|_3.$$

We integrate in time and use the previously established bounds on $\varepsilon^{1/2} \nabla p_\varepsilon$, $\mathrm{div} \, \mathbf{u}_\varepsilon$ and $\mathbf{S} \mathbf{u}_\varepsilon$ (the latter two relying on the bound for $\mathbf{D} \mathbf{u}_\varepsilon$). It remains to estimate $\|V_\varepsilon(p_\varepsilon)\|_3$. To do this, we take $\psi = |V_\varepsilon(p_\varepsilon)|V_\varepsilon(p_\varepsilon)$ in (23) and we obtain

$$\frac{d}{dt} \int_{\Omega} \varepsilon W_{\varepsilon}(p_{\varepsilon}) + \int_{\Omega} \varepsilon W_{\varepsilon}''(p_{\varepsilon}) |\nabla p_{\varepsilon}|^2 + \|V_{\varepsilon}(p_{\varepsilon})\|_{3}^{3} \le \int_{\Omega} 2V_{\varepsilon}(p_{\varepsilon})^2 |\mathbf{S} \mathbf{u}_{\varepsilon}| + \int_{\Omega} V_{\varepsilon}(p_{\varepsilon})^2 |\mathrm{div} \, \mathbf{u}_{\varepsilon}|, \quad (26)$$

where $W_\varepsilon$ is such that $W'_\varepsilon = |V_\varepsilon|V_\varepsilon$ and $W_\varepsilon(0) = 0$. The right-hand side is bounded using the previous estimates. For instance, applying successively Hölder's inequality followed by Young's inequality, we derive

$$\int_{\Omega} V_{\varepsilon}(p_{\varepsilon})^2 |\mathrm{div} \, \mathbf{u}_{\varepsilon}| \le \|V_{\varepsilon}(p_{\varepsilon})\|_{3}^{2} \|\mathrm{div} \, \mathbf{u}_{\varepsilon}\|_{3} \le \frac{1}{3} \|V_{\varepsilon}(p_{\varepsilon})\|_{3}^{3} + \frac{2}{3} \|\mathrm{div} \, \mathbf{u}_{\varepsilon}\|_{3}^{3}.$$

Since $W_\varepsilon$ and $W_\varepsilon''$ are non-negative, and $W_\varepsilon(p_\varepsilon|_{t=0}) = W_\varepsilon(0) = 0$, integration of (26) with respect to time over the interval $(0, T)$ yields the desired bound on $V_\varepsilon(p_\varepsilon)$.

We deduce that there exists a constant $C_2$, independent of $\varepsilon$, such that $\|\varepsilon \partial_t p_\varepsilon\|_{L^2(0,T;W^{-1,2}(\Omega))} \le C_2$. ■

## 3.5 Limit process

By virtue of the bounds obtained in Proposition 3, we deduce the following convergences (up to sub-sequence extractions):
**Proposition 4.** There exists a function $\mathbf{u}$ such that

$$
\begin{align*}
\mathbf{u}_{\varepsilon} \to \mathbf{u} &\quad \text{weakly-\star in } L^{\infty}(0, T; L^2(\Omega)), \\
\nabla \mathbf{u}_{\varepsilon} \to \nabla \mathbf{u} &\quad \text{weakly in } L^3(0, T; L^3(\Omega)), \\
\partial_t \mathbf{u}_{\varepsilon} \to \partial_t \mathbf{u} &\quad \text{weakly in } L^{\frac{3}{2}}(0, T; W^{-1, \frac{3}{2}}(\Omega)), \\
\mathbf{u}_{\varepsilon} \to \mathbf{u} &\quad \text{in } L^3(0, T; L^2(\Omega)).
\end{align*}
$$

There exists a non-negative function $p$ such that

$$
\begin{align*}
p_{\varepsilon} \to p &\quad \text{weakly in } L^{\frac{3}{2}}(0, T; L^{\frac{3}{2}}(\Omega)), \\
\varepsilon \nabla p_{\varepsilon} \to \mathbf{0} &\quad \text{weakly in } L^2(0, T; L^2(\Omega)), \\
\varepsilon \partial_t p_{\varepsilon} \to 0 &\quad \text{weakly in } L^2(0, T; W^{-1,2}(\Omega)).
\end{align*}
$$

There exists a symmetric and traceless tensor function $\sigma$ such that

$$
\sigma_{\varepsilon} \to \sigma \quad \text{weakly in } L^{\frac{3}{2}}(0, T; L^{\frac{3}{2}}(\Omega)).
$$

**Proof -** Weak convergence, up to a sub-sequence, for $\mathbf{u}_\varepsilon$, $\nabla \mathbf{u}_\varepsilon$, $\partial_t \mathbf{u}_\varepsilon$, $p_\varepsilon$ and $\boldsymbol{\sigma}_\varepsilon$ follows from the bounds derived in Proposition 3. Strong convergence on velocity $\mathbf{u}_\varepsilon$ can be deduced from the Aubin-Lions-Simon Lemma (see [11, p.102]).

Since $p_\varepsilon$ weakly converges to $p$, we know that $\partial_t p_\varepsilon$ and $\nabla p_\varepsilon$ converge respectively to $\partial_t p$ and $\nabla p$ in the sense of distributions. Thus, we deduce that $\varepsilon \partial_t p_\varepsilon$ and $\varepsilon \partial_t^2 p_\varepsilon$ converge to zero in the sense of distributions. Since these two sequences are bounded in $L^2(0, T; W^{-1,2}(\Omega))$ and $L^2((0, T) \times \Omega)$, respectively, they converge to zero weakly in these spaces.

The last point that remains to be proven is the non-negativity of the pressure $p$. For $\varphi \in L^3((0, T) \times \Omega)$, $\varphi \ge 0$, we have

$$
\int_0^T \int_{\Omega} p_{\varepsilon} \varphi = \underbrace{\iint_{p_{\varepsilon} \ge 0}}_{A} p_{\varepsilon} \varphi + \underbrace{\iint_{p_{\varepsilon} < 0}}_{B} p_{\varepsilon} \varphi. \qquad (27)
$$

First note that for all $\varepsilon > 0$ we have $A \ge 0$. Moreover, using the Cauchy-Schwarz inequality and the energy estimate (25), we obtain

$$
|B| \le \left( \iint_{p_{\varepsilon} < 0} p_{\varepsilon}^2 \right)^{\frac{1}{2}} \left( \int_0^T \int_{\Omega} \varphi^2 \right)^{\frac{1}{2}} \le \varepsilon^{\frac{1}{2}} \left( \int_0^T \int_{\Omega} p_{\varepsilon} V_{\varepsilon}(p_{\varepsilon}) \right)^{\frac{1}{2}} \left( \int_0^T \int_{\Omega} \varphi^2 \right)^{\frac{1}{2}} \xrightarrow{\varepsilon \to 0} 0.
$$

Finally, since $(p_\varepsilon)_{\varepsilon>0}$ weakly converges to $p$ in $L^{\frac{3}{2}}((0, T) \times \Omega)$, passing to the limit as $\varepsilon \to 0$, equality (27) becomes

$$
\int_0^T \int_{\Omega} p \varphi \ge 0,
$$

which corresponds to the non-negativity of the pressure $p$.

The aim at the end of this section is to demonstrate that we can pass to the limit in relations (22), (23) and (24), especially by handling non-linear terms.

**Step 1: Limit in** (22) Given $\varphi \in L^3(0, T; W_0^{1,3}(\Omega))$, the convergences obtained in Proposition 4 indicate that

$$
\langle \partial_t \mathbf{u}, \varphi \rangle - \int_{\Omega} p \operatorname{div} \varphi + \int_{\Omega} 2|\overline{\mathbf{D}\mathbf{u}}|\overline{|\mathbf{D}\mathbf{u}} : \mathbf{D}\varphi + \int_{\Omega} \boldsymbol{\sigma} : \mathbf{S}\varphi = \langle \mathbf{f}, \varphi \rangle, \quad (28)
$$

where the notation $|\overline{\mathbf{D}\mathbf{u}}|\overline{|\mathbf{D}\mathbf{u}}$ denotes the weak limit of $|\mathbf{D}\mathbf{u}_\varepsilon||\mathbf{D}\mathbf{u}_\varepsilon$ in $L^{\frac{3}{2}}((0, T) \times \Omega)$.
**Step 2: Limit in (23)** In order to perform the limit $\varepsilon \to 0$ in (23), we first note that

$$ \left| \frac{|\mathbf{S}\mathbf{u}_\varepsilon|^2}{|\mathbf{S}\mathbf{u}_\varepsilon| + \varepsilon} - |\mathbf{S}\mathbf{u}_\varepsilon| \right| = \varepsilon \frac{|\mathbf{S}\mathbf{u}_\varepsilon|}{|\mathbf{S}\mathbf{u}_\varepsilon| + \varepsilon} \le \varepsilon \to 0. $$

Thus, for $\psi \in L^{\frac{3}{2}}((0, T) \times \Omega)$ (in practice, this is done for more regular $\psi_k$ functions, typically in $L^2(0, T; W_0^{1,2}(\Omega))$, and which converge to $\psi$), when $\varepsilon$ goes to 0 in (23), we get

$$ \int_{\Omega} \psi \operatorname{div} \mathbf{u} - \int_{\Omega} 2\psi |\overline{\mathbf{S}\mathbf{u}}| + \int_{\Omega} \psi \overline{V(p)} = 0, \quad (29) $$

where $|\overline{\mathbf{S}\mathbf{u}}|$ and $\overline{V(p)}$ respectively correspond to the weak limits of $|\mathbf{S}\mathbf{u}_\varepsilon|$ and $V_\varepsilon(p_\varepsilon)$ in $L^3((0, T) \times \Omega)$.

**Step 3: Limit in (24)** After integrating equation (24) with respect to time we obtain

$$ \frac{1}{2} \|\mathbf{u}_\varepsilon\|_2^2 + \frac{\varepsilon}{2} \|p_\varepsilon\|_2^2 + 4 \int_0^T \|\mathbf{D}\mathbf{u}_\varepsilon\|_3^3 + \varepsilon \int_0^T \|\nabla p_\varepsilon\|_2^2 + \int_0^T \int_\Omega p_\varepsilon V_\varepsilon(p_\varepsilon) = \frac{1}{2} \|\mathbf{u}_{\text{init},\varepsilon}\|_2^2 + \int_0^T \langle \mathbf{f}_\varepsilon, \mathbf{u}_\varepsilon \rangle. $$

In particular, we have

$$ \frac{1}{2} \|\mathbf{u}_\varepsilon\|_2^2 + 4 \int_0^T \|\mathbf{D}\mathbf{u}_\varepsilon\|_3^3 + \int_0^T \int_\Omega p_\varepsilon V_\varepsilon(p_\varepsilon) \le \frac{1}{2} \|\mathbf{u}_{\text{init},\varepsilon}\|_2^2 + \int_0^T \langle \mathbf{f}_\varepsilon, \mathbf{u}_\varepsilon \rangle. \quad (30) $$

Since $(\mathbf{u}_\varepsilon)_{\varepsilon>0}$ converges weakly to $\mathbf{u}$ in $L^\infty(0, T; L^2(\Omega))$, we know that $\liminf_{\varepsilon \to 0} \|\mathbf{u}_\varepsilon\| \ge \|\mathbf{u}\|$. We also use the strong convergence of the sequences $(\mathbf{u}_{\text{init},\varepsilon})_{\varepsilon>0}$ and $(\mathbf{f}_\varepsilon)_{\varepsilon>0}$ to perform the limit $\varepsilon \to 0$ in (30) and obtain

$$ \frac{1}{2} \|\mathbf{u}\|_2^2 + 4 \int_0^T \int_\Omega |\overline{\mathbf{D}\mathbf{u}}|^3 + \int_0^T \int_\Omega \overline{pV(p)} \le \frac{1}{2} \|\mathbf{u}_{\text{init}}\|_2^2 + \int_0^T \langle \mathbf{f}, \mathbf{u} \rangle, \quad (31) $$

where $|\overline{\mathbf{D}\mathbf{u}}|^3$ and $\overline{pV(p)}$ denote the weak limits of $(|\overline{\mathbf{D}\mathbf{u}}|)^3)_{\varepsilon>0}$ and $(p_\varepsilon V_\varepsilon(p_\varepsilon))_{\varepsilon>0}$, respectively, in $L^1((0, T) \times \Omega)$.

## 3.6 Convergence of the non-linear terms and proof of the existence result

In this subsection we will prove that

$$ |\overline{\mathbf{D}\mathbf{u}}| = |\overline{\mathbf{D}\mathbf{u}}|, \quad \overline{V(p)} = \sqrt{p}, \quad |\overline{\mathbf{D}\mathbf{u}}|^3 = |\overline{\mathbf{D}\mathbf{u}}|^3 \quad \text{and} \quad |\overline{\mathbf{S}\mathbf{u}}| = |\overline{\mathbf{S}\mathbf{u}}|. \quad (32) $$

To this end, we make use of the following two lemmas:

### 3.6.1 Convex analysis lemmas

**Lemma 1 (Convexity argument).**

If $H: \mathcal{D} \subset \mathbb{R}^m \to \mathbb{R}$ is a convex and lower semi-continuous function, if $g_\varepsilon: (0, T) \times \Omega \to \mathcal{D}$ weakly converges to $g$ in $L^q((0, T) \times \Omega)$ when $\varepsilon$ tends to 0, if $H(g_\varepsilon): (0, T) \times \Omega \to \mathbb{R}$ weakly converges to $\overline{H(g)}$ in $L^r((0, T) \times \Omega)$ when $\varepsilon$ tends to 0 then we have

$$ H(g) \le \overline{H(g)}. $$

**Proof:** This is a classical result that can be seen as a consequence of the equivalence of the notion of convex closed set for weak and strong topologies, see [13, p. 38]. For any $\varphi \in L^{r'}((0, T) \times \Omega)$, $\varphi \ge 0$ where $\frac{1}{r} + \frac{1}{r'} = 1$, the function

$$ \Phi: h \in L^q((0, T) \times \Omega) \mapsto \int_0^T \int_\Omega H(h) \varphi $$

is convex and lower semi-continuous. Following [13, Corollaire III.8] we deduce that if $g_\varepsilon$ weakly converges to $g$ in $L^q((0, T) \times \Omega)^m$ then $\Phi(g) \le \liminf_{\varepsilon \to 0} \Phi(g_\varepsilon)$. We deduce that for all $\varphi \in L^{r'}((0, T) \times \Omega)$, $\varphi \ge 0$, we have

$$ \int_0^T \int_\Omega H(g) \varphi \le \int_0^T \int_\Omega \overline{H(g)} \varphi, $$
which corresponds to the announced result.

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

Clearly, this result can also be extended to the case where the functions $H$ is concave and upper semi-continuous. In this case, the function $-H$ satisfies the assumptions of Lemma 1 and we obtain

$$H(g) \geq \overline{H(g)}.$$

**Corollary 1.** We have the following inequalities

$$\overline{pV(p)} \geq p\sqrt{p} \quad \text{and} \quad \overline{V(p)} \leq \sqrt{p}.$$

**Proof of the first inequality:** Let $\varphi \in L^{\infty}((0, T) \times \Omega)$ with $\varphi \geq 0$. For all $\varepsilon > 0$ we have

$$\int_{0}^{T} \int_{\Omega} p_{\varepsilon} V_{\varepsilon}(p_{\varepsilon}) \varphi = \underbrace{\iint_{p_{\varepsilon} \leq 0} p_{\varepsilon} V_{\varepsilon}(p_{\varepsilon}) \varphi}_{A} + \underbrace{\iint_{p_{\varepsilon} > 0} p_{\varepsilon}(V_{\varepsilon}(p_{\varepsilon}) - \sqrt{p_{\varepsilon}}) \varphi}_{B} + \underbrace{\iint_{p_{\varepsilon} > 0} p_{\varepsilon} \sqrt{p_{\varepsilon}} \varphi}_{C}.$$

- Since for all $x \leq 0$ we have $xV_{\varepsilon}(x) \geq 0$ we deduce that $A \geq 0$.

- Since for all $x > 0$ we have $|V_{\varepsilon}(x) - \sqrt{x}| \leq \varepsilon$, we deduce that $\lim_{\varepsilon \to 0} B = 0$.

- Finally, applying Lemma 1 with the convex and continuous function $H: x \in \mathbb{R}_{+} \mapsto x\sqrt{x}$, and using the fact that the limit pressure $p$ is non-negative, we deduce that $\liminf_{\varepsilon \to 0} C \geq \int_{0}^{T} \int_{\Omega} p\sqrt{p} \varphi$.

We then conclude that

$$\int_{0}^{T} \int_{\Omega} \overline{pV(p)} \varphi \geq \int_{0}^{T} \int_{\Omega} p\sqrt{p} \varphi.$$

**Proof of the second inequality:** Let $\varphi \in L^{\frac{3}{2}}((0, T) \times \Omega)$ with $\varphi \geq 0$. For all $\varepsilon > 0$ we have

$$\int_{0}^{T} \int_{\Omega} V_{\varepsilon}(p_{\varepsilon}) \varphi = \underbrace{\iint_{p_{\varepsilon} \leq 0} V_{\varepsilon}(p_{\varepsilon}) \varphi}_{A} + \underbrace{\iint_{p_{\varepsilon} > 0} V_{\varepsilon}(p_{\varepsilon}) \varphi}_{B}.$$

- Since for all $x \leq 0$ we have $V_{\varepsilon}(x) \leq 0$, we deduce $A \leq 0$.

- Since for all $x > 0$ we have $V_{\varepsilon}(x) \leq \sqrt{x}$, we deduce $B \leq \iint_{p_{\varepsilon} > 0} \sqrt{p_{\varepsilon}} \varphi$. Then, applying Lemma 1 (in its concave version, as detailed immediately following the lemma's proof) with the concave and continuous function $H: x \in \mathbb{R}_{+} \mapsto \sqrt{x}$, and using the fact that the limit pressure $p$ is non-negative, we deduce that $\liminf_{\varepsilon \to 0} B \leq \int_{0}^{T} \int_{\Omega} \sqrt{p} \varphi$.

We then conclude that

$$\int_{0}^{T} \int_{\Omega} \overline{V(p)} \varphi \leq \int_{0}^{T} \int_{\Omega} \sqrt{p} \varphi,$$

that completes the proof of Corollary 1.

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

The following lemma can be found in [29, Lemma 1], see also [17, p.1715 and p.1730]:

**Lemma 2.** If $Du_{\varepsilon} \to Du$ weakly in $L^3((0, T) \times \Omega)$ and $|Du_{\varepsilon}|Du_{\varepsilon} \to \overline{|Du|Du}$ weakly in $L^{\frac{3}{2}}((0, T) \times \Omega)$ then

$$\int_{0}^{T} \int_{\Omega} \left( \overline{|Du|Du} : Du - 2|Du|^3 \right) \leq 0.$$

where $|Du|^3$ refers to the weak limit of $|Du_{\varepsilon}|^3$ in $L^1((0, T) \times \Omega)$.

Furthermore, if there is equality, then $|Du|Du = |Du|Du$.
### 3.6.2 Existence proof

With the help of Lemmas 1 and 2, we can conclude the existence proof stated in the theorem by choosing $\varphi = \mathbf{u}$ in equation (28) and $\psi = p$ in equation (29), integrating both equations over time, and summing the resulting expressions. This yields

$$ \frac{1}{2} \|\mathbf{u}\|_2^2 + \int_0^T \int_{\Omega} 2|\overline{\mathbf{Du}}|\overline{\mathbf{Du}} : \mathbf{Du} + \int_0^T \int_{\Omega} p\overline{V(p)} + \int_0^T \int_{\Omega} (\boldsymbol{\sigma} : \mathbf{Su} - 2p|\overline{\mathbf{Su}}|) = \frac{1}{2} \|\mathbf{u}_{\text{init}}\|_2^2 + \int_0^T \langle \mathbf{f}, \mathbf{u} \rangle. $$

By subtracting equation (31), we obtain

$$ \underbrace{\int_0^T \int_{\Omega} 2(\overline{|\mathbf{Du}|\mathbf{Du}} : \mathbf{Du} - 2|\overline{\mathbf{Du}}|^3)}_{A} + \underbrace{\int_0^T \int_{\Omega} (p\overline{V(p)} - \overline{pV(p)})}_{B} + \underbrace{\int_0^T \int_{\Omega} (\boldsymbol{\sigma} : \mathbf{Su} - 2p|\overline{\mathbf{Su}}|)}_{C} \geq 0. \quad (33) $$

- By Lemma 2, we have $A \leq 0$.

- From Corollary 1, we directly obtain $B \leq 0$.

- By the convexity of $x \mapsto |x|$ and Lemma 1, we have $|\mathbf{Su}| \leq |\overline{\mathbf{Su}}|$. Moreover, we have $|\boldsymbol{\sigma}| \leq p$. Consequently, by applying the Cauchy-Schwarz inequality we write

$$ \boldsymbol{\sigma} : \mathbf{Su} \leq 2|\boldsymbol{\sigma}| |\overline{\mathbf{Su}}| \leq 2p|\overline{\mathbf{Su}}| \leq 2p|\overline{\mathbf{Su}}|. \quad (34) $$

which implies that $C \leq 0$.

Since $A + B + C \geq 0$ with $A \leq 0$, $B \leq 0$ and $C \leq 0$, it follows that $A = B = C = 0$.

- According to the conclusion of Lemma 2, the equality $A = 0$ implies $|\overline{\mathbf{Du}}|\overline{\mathbf{Du}} = |\mathbf{Du}|\mathbf{Du}$.

- Similarly, the equality $B = 0$ implies $\overline{V(p)} = \sqrt{p}$.

- Finally, even if this point is not useful in the proof of Theorem 1, we can note that $C = 0$ implies $\boldsymbol{\sigma} : \mathbf{Su} = 2p|\overline{\mathbf{Su}}|$ so that equality holds in (34), implying that $\boldsymbol{\sigma} : \mathbf{Su} = 2p|\overline{\mathbf{Su}}|$.

In particular, this last point proves that $2p|\overline{\mathbf{Su}}| = 2p|\overline{\mathbf{Su}}|$ but does not prove that $|\overline{\mathbf{Su}}| = |\overline{\mathbf{Su}}|$. To obtain the latter result, we proceed in the following four steps.

**Step 1:** $|\overline{\mathbf{Du}}|^3 = |\mathbf{Du}|^3$

This is a direct consequence of Lemma 2. Indeed, since $A = 0$ and $|\overline{\mathbf{Du}}|\overline{\mathbf{Du}} = |\mathbf{Du}|\mathbf{Du}$, we have

$$ \int_0^T \int_{\Omega} (|\mathbf{Du}|^3 - |\overline{\mathbf{Du}}|^3) = 0. $$

Moreover the function $x \mapsto x^3$ being convex, we apply Lemma 1 to obtain $|\overline{\mathbf{Du}}|^3 - |\mathbf{Du}|^3 \leq 0$. We deduce the equality $|\overline{\mathbf{Du}}|^3 = |\mathbf{Du}|^3$.

**Step 2:** $|\overline{\mathbf{Du}}|^2 = |\mathbf{Du}|^2$

Using the convexity of $x \mapsto x^2$ and the concavity of $x \mapsto x^{\frac{2}{3}}$, we deduce from Lemma 1 that

$$ |\overline{\mathbf{Du}}|^2 \leq |\overline{\mathbf{Du}}|^2 = (|\overline{\mathbf{Du}}|^3)^{\frac{2}{3}} \leq |\overline{\mathbf{Du}}|^3 \stackrel{\frac{2}{3}}{=} |\mathbf{Du}|^2. $$

All these terms are therefore equal, and in particular $|\overline{\mathbf{Du}}|^2 = |\mathbf{Du}|^2$.

**Step 3:** $|\overline{\mathbf{Su}}|^2 = |\mathbf{Su}|^2$

We use the identity

$$ |\overline{\mathbf{Du}}_{\varepsilon}|^2 = |\mathbf{Su}_{\varepsilon}|^2 + \frac{1}{6}|\operatorname{div} \mathbf{u}_{\varepsilon}|^2. \quad (35) $$

Let $\xi \geq 0$ and $\zeta \geq 0$ be the elements of $L^2((0, T) \times \Omega)$ such that $|\overline{\mathbf{Su}}|^2 = |\mathbf{Su}|^2 + \xi$ and $|\operatorname{div} \mathbf{u}|^2 = |\operatorname{div} \mathbf{u}|^2 + \zeta$ then we pass to the limit in (35) using the fact that $|\overline{\mathbf{Du}}|^2 = |\mathbf{Du}|^2$:

$$ |\overline{\mathbf{Du}}|^2 = |\mathbf{Su}|^2 + \xi + \frac{1}{6}|\operatorname{div} \mathbf{u}|^2 + \frac{1}{6}\zeta. $$
Since we also have $|\mathbf{Du}|^2 = |\mathbf{Su}|^2 + \frac{1}{6}|\mathrm{div} \, \mathbf{u}|^2$, it follows that $\xi + \frac{1}{6}\zeta = 0$. Given the positivity of $\xi$ and $\zeta$, we conclude that $\xi = \zeta = 0$, and in particular $|\overline{\mathbf{Su}}|^2 = |\mathbf{Su}|^2$.

**Step 4:** $|\overline{\mathbf{Su}}| = |\mathbf{Su}|$

Using the convexity of $x \mapsto |x|$, and the concavity of $x \mapsto x^{\frac{1}{2}}$, we deduce from Lemma 1 that

$$|\mathbf{Su}| \leq |\overline{\mathbf{Su}}| = (|\mathbf{Su}|^2)^{\frac{1}{2}} \leq |\overline{\mathbf{Su}}|^{\frac{2}{2}} = |\mathbf{Su}|.$$

All these terms are therefore equal, and in particular $|\overline{\mathbf{Su}}| = |\mathbf{Su}|$.

To conclude the proof of Theorem 1, simply note that equations (28), (29) and (31) exactly correspond to Definition 1 of a weak solution since we have proved (32). ■

## 3.7 Proof of the uniqueness result

The proof of the uniqueness of a weak solution is classical. Consider two weak solutions $(\mathbf{u}_1, p_1, \sigma_1)$ and $(\mathbf{u}_2, p_2, \sigma_2)$. We write the equations (13) and (14) for both solutions, and then take the difference. Denoting $\mathbf{u} = \mathbf{u}_2 - \mathbf{u}_1$, $p = p_2 - p_1$ and $\sigma = \sigma_2 - \sigma_1$, we obtain, for all $\varphi \in L^3(0, T; W_0^{1,3}(\Omega))$ and $\psi \in L^{\frac{3}{2}}((0, T) \times \Omega)$:

$$\langle \partial_t \mathbf{u}, \varphi \rangle - \int_{\Omega} p \mathrm{div} \, \varphi + \int_{\Omega} 2(|\mathbf{Du}_2|\mathbf{Du}_2 - |\mathbf{Du}_1|\mathbf{Du}_1) : \mathbf{D}\varphi + \int_{\Omega} \sigma : \mathbf{S}\varphi = 0, \\ \int_{\Omega} \psi \mathrm{div} \, \mathbf{u} - \int_{\Omega} 2\psi(|\mathbf{Su}_2| - |\mathbf{Su}_1|) + \int_{\Omega} \psi(\sqrt{p_2} - \sqrt{p_1}) = 0.$$

We then choose $\varphi = \mathbf{u}$ and $\psi = p$, and add the resulting equations:

$$\frac{1}{2} \frac{d}{dt} \int_{\Omega} |\mathbf{u}|^2 + \int_{\Omega} \underbrace{(|\mathbf{Du}_2|\mathbf{Du}_2 - |\mathbf{Du}_1|\mathbf{Du}_1) : (\mathbf{Du}_2 - \mathbf{Du}_1)}_{A} + \int_{\Omega} \underbrace{(\sqrt{p_2} - \sqrt{p_1})(p_2 - p_1)}_{B} + \int_{\Omega} \underbrace{(\sigma_2 - \sigma_1) : (\mathbf{Su}_2 - \mathbf{Su}_1) - 2(p_2 - p_1)(|\mathbf{Su}_2| - |\mathbf{Su}_1|)}_{C} \leq 0. \quad (36)$$

We remark that $A \geq (|\mathbf{Du}_2| + |\mathbf{Du}_1|)(|\mathbf{Du}_2| - |\mathbf{Du}_1|)^2 \geq 0$ and that $B = (\sqrt{p_2} + \sqrt{p_1})(\sqrt{p_2} - \sqrt{p_1})^2 \geq 0$. Moreover, using Proposition 1 and the Cauchy-Schwarz inequality, we also have $C \geq 0$. Consequently, estimate (36) implies that $\|\mathbf{u}\|_2$ is a positive and non-increasing function of time. Since $\mathbf{u}_2$ and $\mathbf{u}_1$ satisfy the same initial condition, we deduce that $\mathbf{u} = 0$.

We next deduce from the estimate (36) that $A = B = C = 0$. In particular, the equality $B = 0$ implies that $p = 0$. ■

# 4 Adaptation to more general models

## 4.1 Taking into account the variation in volume fraction

As explained in the introduction, the complete granular flow model also describes the evolution of the grain volume fraction, denoted $\phi$. Although the previous results have overcome a number of difficulties relating to complex flows (non-linear rheology and dilatation phenomena), the complete model given by (1)-(4) in the introduction presents additional difficulties. The key point is to specify the functions $\alpha$ and $\beta$ which depend on the volume fraction $\phi$. More specifically, it is important to understand how to describe the expansion law expressing $\mathrm{div} \, \mathbf{u}$. From Roux-Radjaï's work [27], we know that this law can be written as

$$\mathrm{div} \, \mathbf{u} = 2K|\mathbf{Su}|(I - I_{\mathrm{eq}}(\phi)) \quad \text{where} \quad I = \frac{2d|\mathbf{Su}|}{\sqrt{p/\rho_0}}.$$
The non-dimensional number I denotes the inertial number (depending on the diameter d of the grains, and on their density ρ₀) and Ieq(φ) its equilibrium value. For the latter, several examples are given in the literature, see [12, 19, 21, 25], and we choose [28, p.929]:

$$I_{\text{eq}}(\phi) = \frac{\phi_{\text{max}} - \phi}{\phi - \phi_{\text{min}}}$$

where the constants 0 < φmin < φmax < 1 corresponds to the "extreme" values of the volume fraction φ. This choice will ensure mathematically that φ remains effectively bounded between φmin and φmax, see Proposition 6 below. In order to ensure stability of the model (see [4, 5, 6, 14]), we choose the proportionality coefficient as K = φ-φmin/δφI with δφ = φmax - φmin, so that the dilatation law is written as

$$\operatorname{div} \mathbf{u} = 2 \frac{\phi - \phi_{\min}}{\delta \phi} |\mathbf{S}\mathbf{u}| - \frac{\phi_{\max} - \phi}{\delta \phi} \frac{\sqrt{p}}{d \sqrt{\rho_0}}.$$

Regarding the mathematical analysis, we propose here a first theoretical result on the existence and uniqueness for a model that is relatively close to the full model. This considered model is as follows:

$$\mathcal{D}_t(\phi, \mathbf{u}) + \nabla p - \operatorname{div}(2|\mathbf{D}\mathbf{u}|\mathbf{D}\mathbf{u}) = \mathbf{f} + \operatorname{div} \boldsymbol{\sigma}, \qquad (37)$$

$$\boldsymbol{\sigma} : \mathbf{S}\mathbf{u} = 2(\phi - \phi_{\min})p|\mathbf{S}\mathbf{u}|, \quad |\boldsymbol{\sigma}| \le (\phi - \phi_{\min})p, \quad \boldsymbol{\sigma} = \boldsymbol{\sigma}^{\top} \quad \text{and} \quad \operatorname{tr} \boldsymbol{\sigma} = 0, \qquad (38)$$

$$\operatorname{div} \mathbf{u} = 2(\phi - \phi_{\min})|\mathbf{S}\mathbf{u}| - (\phi_{\max} - \phi)\sqrt{p}, \qquad (39)$$

$$\partial_t \phi + \operatorname{div}(\phi \mathbf{u}) = \xi \mathcal{H}, \qquad (40)$$

where

$$\mathcal{D}_t(\phi, \mathbf{u}) = \frac{1}{2} \left( (\partial_t(\phi \mathbf{u}) + \operatorname{div}(\phi \mathbf{u} \otimes \mathbf{u})) + (\phi(\partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u})) \right).$$

In practice, the only difference with the full model lies in the additional term ξH where ξ is a real parameter, fixed and positive (possibly very small). This term H will be chosen in the following form

$$\mathcal{H} = \Delta\phi - \phi\sqrt{p}.$$

The contribution of Δφ ensures compactness for sequence (φε)ε>0 satisfying such a problem, and gives sense to the non-linear terms such as (φmax - φ)√p. Furthermore, the contribution of φ√p ensures that φ ≤ φmax - ξ, thereby allowing the energy estimate to yield a bound on the pressure.

**Remark 5.** The form of the time derivative as the average between the convected derivative and its conservative form ensures an energy estimate independent of the evolution equation (40) on φ. In the case where ξ = 0, both forms are equivalent. At this stage, we do not know how to obtain sufficient regularity for the solutions to assert the existence of a solution when ξ = 0. In particular, we do not know whether the solutions constructed here converge to a solution when ξ goes to 0.

**Theorem 2.** Let ξ > 0, uinit ∈ L²(Ω), φinit ∈ L²(Ω) and f ∈ L³/2(0, T; W⁻¹,³/2(Ω)). If φmin ≤ φinit ≤ φmax - ξ then there exists a weak solution of (37)–(40) satisfying the initial conditions u|t=0 = uinit and φ|t=0 = φinit.

The proof follows the same arguments as in the case of Theorem 1. The only difference lies in the control of φ. We must therefore ensure that there is always an energy estimate that, in addition to providing control over u and p, also allows us to control φ. To avoid repeating the full proof of the previous theorem, we simply derive the energy estimate in a formal manner.

**Proposition 5.** Any regular solution (u, p, σ, φ) to (37)–(40) satisfies

$$\frac{d}{dt} \int_{\Omega} \phi \frac{|\mathbf{u}|^2}{2} + 4 \int_{\Omega} |\mathbf{D}\mathbf{u}|^3 + \int_{\Omega} (\phi_{\max} - \phi) p \sqrt{p} = \int_{\Omega} \mathbf{f} \cdot \mathbf{u}. \qquad (41)$$

**Proof** - The result essentially follows by taking the scalar product of equation (37) with u. We note that, independently of equation (40), we have

$$\mathcal{D}_t(\phi, \mathbf{u}) \cdot \mathbf{u} = \partial_t \left( \phi \frac{|\mathbf{u}|^2}{2} \right) + \operatorname{div} \left( \phi \frac{|\mathbf{u}|^2}{2} \mathbf{u} \right).$$
Equation (39) allows us to control the pressure term

$$\nabla p \cdot \mathbf{u} = \operatorname{div}(p\mathbf{u}) - p\operatorname{div}\mathbf{u} = \operatorname{div}(p\mathbf{u}) - 2(\phi - \phi_{\min})p|\mathbf{u}|^2 + (\phi_{\max} - \phi)p\sqrt{p},$$

while equation (38) is useful for cancelling out the stress part:

$$\operatorname{div} \boldsymbol{\sigma} \cdot \mathbf{u} = \operatorname{div}(\boldsymbol{\sigma} \cdot \mathbf{u}) - \boldsymbol{\sigma} : \mathbf{S} \mathbf{u} = \operatorname{div}(\boldsymbol{\sigma} \cdot \mathbf{u}) - 2(\phi - \phi_{\min})p|\mathbf{S} \mathbf{u}|.$$

After integration over $\Omega$, the sum of these contributions provides the result stated in Proposition 5.

To use equation (41), we need to prove that the terms containing $\phi|\mathbf{u}|^2$ and $(\phi_{\max} - \phi)p\sqrt{p}$ are positive, more precisely that $0 < \phi \le \phi_{\max}$. In addition, if we want to obtain an estimate on the pressure, we need to ensure that $\phi_{\max} - \phi$ does not vanish.

**Proposition 6.** Any regular solution $(\mathbf{u}, p, \boldsymbol{\sigma}, \phi)$ to (37)-(40) such that $\phi_{\min} \le \phi|_{t=0} \le \phi_{\max} - \xi$ satisfies

$$\phi_{\min} \le \phi \le \phi_{\max} - \xi, \qquad (42)$$

and

$$\frac{d}{dt} \|\phi\|_2^2 + 2\xi \|\nabla \phi\|_2^2 \le \phi_{\max}^2 \int_{\Omega} (|\operatorname{div} \mathbf{u}| + 2\sqrt{p}). \qquad (43)$$

**Proof – Step 1: Upper bound** First note that equation (40) implies that, for $\beta: \mathbb{R} \to \mathbb{R}$, we have

$$\partial_t \beta(\phi) + \operatorname{div}(\beta(\phi)\mathbf{u}) + (\phi\beta'(\phi) - \beta(\phi))\operatorname{div}\mathbf{u} = \xi\beta'(\phi)\Delta\phi - \xi\phi\beta'(\phi)\sqrt{p}, \qquad (44)$$

In order to achieve the upper bound for $\phi$, we choose the function $\beta$ defined by

$$\beta(\phi) = \begin{cases} 0 & \text{if } \phi < \phi_{\max} - \xi, \\ \phi - (\phi_{\max} - \xi) & \text{if } \phi \ge \phi_{\max} - \xi. \end{cases}$$

Integrating (44) over $\Omega$, and thanks to this choice of the function $\beta$, we obtain

$$\frac{d}{dt} \int_{\Omega} \beta(\phi) + \int_{\mathcal{E}_+} (\phi_{\max} - \xi) \operatorname{div} \mathbf{u} = \xi \int_{\mathcal{E}_+} \Delta\phi - \int_{\mathcal{E}_+} \xi\phi\sqrt{p}, \qquad (45)$$

where $\mathcal{E}_+ = \{x \in \Omega; \phi(t, x) \ge \phi_{\max} - \xi\}$. The unit outgoing normal vector at $\mathcal{E}_+$ is given by $-\nabla\phi/|\nabla\phi|$ so that Stokes' formula allows us to write

$$\xi \int_{\mathcal{E}_+} \Delta \phi = -\xi \int_{\partial \mathcal{E}_+} |\nabla \phi| \le 0.$$

Moreover, the expression for the divergence of the velocity given by equation (39) indicates that

$$\phi \ge \phi_{\max} - \xi \implies \operatorname{div} \mathbf{u} + \xi\sqrt{p} \ge 0$$

We then deduce that

$$\int_{\mathcal{E}_+} (\phi_{\max} - \xi) \operatorname{div} \mathbf{u} + \int_{\mathcal{E}_+} \xi\phi\sqrt{p} \ge \int_{\mathcal{E}_+} (\phi_{\max} - \xi)(\operatorname{div} \mathbf{u} + \xi\sqrt{p}) \ge 0.$$

Consequently, equation (45) implies $\frac{d}{dt} \int \beta(\phi) \le 0$. Thus, if $\phi|_{t=0} \le \phi_{\max} - \xi$, i.e. $\beta(\phi|_{t=0}) = 0$, then $\beta(\phi) = 0$ so the result follows, namely $\phi \le \phi_{\max} - \xi$.

**Step 2: Lower bound** In the same spirit, to obtain the lower bound, we use equation (44) together with the function $\beta$, defined by

$$\beta(\phi) = \begin{cases} 0 & \text{if } \phi > \phi_{\min}, \\ \phi_{\min} - \phi & \text{if } \phi \le \phi_{\min}. \end{cases}$$

By integrating (44) over $\Omega$, we get

$$\frac{d}{dt} \int_{\Omega} \beta(\phi) - \int_{\mathcal{E}_-} \phi_{\min} \operatorname{div} \mathbf{u} = -\xi \int_{\mathcal{E}_-} \Delta\phi + \int_{\mathcal{E}_-} \xi\phi\sqrt{p}, \qquad (46)$$
where $\mathcal{E}_- = \{x \in \Omega; \phi(t, x) \le \phi_{\min}\}$. This time, the unit outgoing normal vector at $\mathcal{E}_-$ is given by $\nabla\phi/|\nabla\phi|$ so that Stokes' formula allows us to write

$$-\xi \int_{\mathcal{E}_-} \Delta\phi = -\xi \int_{\partial\mathcal{E}_-} |\nabla\phi| \le 0.$$

The expression for the divergence of the velocity in equation (39) gives

$$\phi \le \phi_{\min} \quad \Rightarrow \quad \operatorname{div} \mathbf{u} \le -(\phi_{\max} - \phi_{\min})\sqrt{p}.$$

Since $\xi \le \phi_{\max} - \phi_{\min}$, we obtain the following inequality

$$\int_{\mathcal{E}_-} \xi \phi \sqrt{p} + \int_{\mathcal{E}_-} \phi_{\min} \operatorname{div} \mathbf{u} \le \int_{\mathcal{E}_-} \phi_{\min} (\xi - (\phi_{\max} - \phi_{\min})) \sqrt{p} \le 0.$$

Consequently, equation (46) implies that $\frac{d}{dt} \int \beta(\phi) \le 0$. Thus, if we have $\phi|_{t=0} \ge \phi_{\min}$, i.e. $\beta(\phi|_{t=0}) = 0$, then $\beta(\phi) = 0$ so that the result follows, namely $\phi \ge \phi_{\min}$.

**Step 3: $H^1$ estimate** By multiplying equation (40) by $2\phi$ and integrating over $\Omega$, we obtain

$$\frac{d}{dt} \int_{\Omega} |\phi|^2 + 2\xi \int_{\Omega} |\nabla\phi|^2 = - \int_{\Omega} 2\phi \operatorname{div}(\phi\mathbf{u}) - \int_{\Omega} 2\phi^2 \sqrt{p}.$$

Using integration by parts, we can write the right-hand side in the form

$$\frac{d}{dt} \|\phi\|_2^2 + 2\xi \|\nabla\phi\|_2^2 = - \int_{\Omega} \phi^2 (\operatorname{div} \mathbf{u} + 2\sqrt{p}).$$

The bound $|\phi| \le \phi_{\max}$ allows to conclude.

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

**Ideas for the proof of Theorem 2** The method is the same as the one described in the proof of Theorem 1. We construct a sequence of solutions $(\mathbf{u}_\epsilon, p_\epsilon, \sigma_\epsilon, \phi_\epsilon)_{\epsilon>0}$ of an approximate problem, and we show that this solution converges to the solution of the problem (37)-(40) when $\epsilon$ tends to 0.

We focus on estimating the volume fraction $\phi_\epsilon$, since the other quantities are bounded similarly to those in the proof of Theorem 1. The key point is to ensure that the sequence $(\phi_\epsilon)_{\epsilon>0}$ converges strongly to a solution $\phi$, which will enable us to pass to the limit in all terms.

Given that $\phi_{\max} - \phi_\epsilon \ge \xi > 0$, estimate (41) ensures that the sequences $(\operatorname{div} \mathbf{u}_\epsilon)_{\epsilon>0}$ and $(\sqrt{p_\epsilon})_{\epsilon>0}$ are uniformly bounded in $L^3((0, T) \times \Omega)$ with respect to $\epsilon$. From estimate (43), we deduce that $(\phi_\epsilon)_{\epsilon>0}$ is uniformly bounded in $L^2(0, T; H^1(\Omega))$.

To obtain compactness, we consider the sequence $(\partial_t \phi_\epsilon)_{\epsilon>0}$. By writing

$$\partial_t \phi = \xi \Delta \phi - \operatorname{div}(\phi \mathbf{u}) - \phi \sqrt{p},$$

we deduce from the previous estimates that $(\partial_t \phi_\epsilon)_{\epsilon>0}$ is uniformly bounded in $L^2(0, T; H^{-1}(\Omega))$. Finally, the Aubin-Lions-Simon Lemma (see [11, p.102]) leads to the conclusion that $(\phi_\epsilon)_{\epsilon>0}$ converges strongly to $\phi$ in $L^2((0, T) \times \Omega)$. ■

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

**Remark 6.** Propositions 5 and 6 are true even if $\xi = 0$. However, the proof of Theorem 2 is no longer correct. Indeed, when $\xi = 0$, we don't have $H^1$ estimate for $\phi$ (see relation (43)), nor an estimate for $p$ (see relation (41), especially when $\phi_{\max} - \phi$ vanishes). To prove such a theorem when $\xi = 0$, we therefore need to find other arguments.

## 4.2 $\mu(I)$-rheology

The $\mu(I)$-rheology relates stress, pressure and shear in a manner similar to equation (5). The main difference is that the stress threshold, corresponding to $\alpha(\phi)p$ in (5), is written $\mu(I)p$, the function $\mu$ being experimentally given with respect to the inertial number $I$, see for example [3]. In this context, if Barker's stability
conditions [28] are to be satisfied, the dilation relation is also imposed, and involves I-dependent functions. By following these principles the following model is proposed in [14]:

$$ \phi\rho_0(\partial_t\mathbf{u} + \mathbf{u} \cdot \nabla\mathbf{u}) + \nabla p = \phi\rho_0\mathbf{g} + \operatorname{div}\boldsymbol{\sigma}, \quad (47) $$

$$ \boldsymbol{\sigma} : \mathbf{S}\mathbf{u} = 2\mu(I)p|\mathbf{S}\mathbf{u}|, \quad |\boldsymbol{\sigma}| \le \mu(I)p, \quad \boldsymbol{\sigma} = \boldsymbol{\sigma}^{\mathrm{T}} \quad \text{and} \quad \operatorname{tr}\boldsymbol{\sigma} = 0, \quad (48) $$

$$ \operatorname{div}\mathbf{u} = 2F(I)|\mathbf{S}\mathbf{u}| - \gamma I_{\mathrm{eq}}(\phi)F(I_{\mathrm{eq}}(\phi))\sqrt{p}, \quad (49) $$

$$ \partial_t\phi + \operatorname{div}(\phi\mathbf{u}) = 0, \quad (50) $$

where $\gamma = 2/d\sqrt{\rho_0}$ and where the dimensionless functions $\mu$ and $F$ are explicitly given (see [14] for more details). In particular, these functions are smooth, continuous, and equal at 0. We denote $\alpha_0 = \mu(0) = F(0)$.

**Dimensionless procedure** In order to rewrite the system in dimensionless form, we apply the following change of variables and unknowns:

$$ t = T\tilde{t}, \quad x = L\tilde{x}, \quad \mathbf{u} = U\tilde{\mathbf{u}}, \quad p = gL\rho_0\tilde{p}, \quad \boldsymbol{\sigma} = gL\rho_0\tilde{\boldsymbol{\sigma}}, $$

where $T$, $L$ and $U$ represent the characteristic time, length and velocity scales, respectively. Introducing the dimensionless numbers $\varepsilon$, $\mathfrak{F}\mathbf{r}$ and $\mathfrak{D}\mathbf{i}$ defined as follows

$$ \varepsilon = \frac{TU}{L}, \quad \mathfrak{F}\mathbf{r}^2 = \frac{U^2}{gL}, \quad \mathfrak{D}\mathbf{i} = \frac{d}{L}\mathfrak{F}\mathbf{r}, $$

yields the following system where the tilde notation has been omitted for conciseness:

$$ \phi(\partial_t \mathbf{u} + \varepsilon \mathbf{u} \cdot \nabla \mathbf{u}) + \frac{\varepsilon}{\mathfrak{F}\mathbf{r}^2} \nabla p = - \frac{\varepsilon}{\mathfrak{F}\mathbf{r}^2} \phi \mathbf{e} + \frac{\varepsilon}{\mathfrak{F}\mathbf{r}^2} \operatorname{div} \boldsymbol{\sigma}, \quad (51) $$

$$ \boldsymbol{\sigma} : \mathbf{S}\mathbf{u} = 2\mu(\mathfrak{D}\mathbf{i}I)p|\mathbf{S}\mathbf{u}|, \quad |\boldsymbol{\sigma}| \le \mu(\mathfrak{D}\mathbf{i}I)p, \quad \boldsymbol{\sigma} = \boldsymbol{\sigma}^{\mathrm{T}} \quad \text{and} \quad \operatorname{tr}\boldsymbol{\sigma} = 0, \quad (52) $$

$$ \operatorname{div}\mathbf{u} = 2F(\mathfrak{D}\mathbf{i}I)|\mathbf{S}\mathbf{u}| - \frac{1}{\mathfrak{D}\mathbf{i}}I_{\mathrm{eq}}(\phi)F(I_{\mathrm{eq}}(\phi))\sqrt{p}, \quad (53) $$

$$ \partial_t\phi + \varepsilon\operatorname{div}(\phi\mathbf{u}) = 0. \quad (54) $$

Note that in this dimensionless form, the inertial number $I$ is written as: $I = \frac{2|\mathbf{S}\mathbf{u}|}{\sqrt{p}}$.

**Asymptotic problem** We aim to study this system when the volume fraction $\phi$ is nearly constant and close to the maximum packing fraction $\phi_{\max}$. Choosing $\varepsilon$ as a small parameter, let us introduce

$$ \phi = \phi_0 + \varepsilon\psi, $$

where $\phi_0$ is a constant close to $\phi_{\max}$ in the following sense: $\phi_{\max} - \phi_0 = \mathcal{O}(\varepsilon)$. More generally, we write

$$ \mathbf{u} = \mathbf{v} + \mathcal{O}(\varepsilon), \quad p = q + \mathcal{O}(\varepsilon) \quad \text{and} \quad \boldsymbol{\sigma} = \boldsymbol{\tau} + \mathcal{O}(\varepsilon). $$

To retain both rheological and expansion-related contributions, we also assume that

$$ \mathfrak{F}\mathbf{r}^2 = \mathcal{O}(\varepsilon) \quad \text{and} \quad \mathfrak{D}\mathbf{i} = \mathcal{O}(\varepsilon). $$

A typical illustration of this situation is shown below:

$$ L = 10^{-1} \, \mathrm{m}, \quad U = 10^{-1} \, \mathrm{m.s}^{-1}, \quad T = 10^{-2} \, \mathrm{s}, \quad d = 10^{-2} \, \mathrm{m} \quad \text{and} \quad g = 10 \, \mathrm{m.s}^{-2}. $$

In this case, we have $\varepsilon = 10^{-2}$, $\mathfrak{F}\mathbf{r}^2 = 10^{-2}$ and $\mathfrak{D}\mathbf{i} = 10^{-2}$. By keeping only the terms of leading order in $\varepsilon$ in the system (51)–(54), we obtain

$$ \phi_0 \partial_t \mathbf{v} + \lambda_0 \nabla q = -\phi_0 \lambda_0 \mathbf{e} + \lambda_0 \operatorname{div} \boldsymbol{\tau}, \quad (55) $$

$$ \boldsymbol{\tau} : \mathbf{S}\mathbf{v} = 2\alpha_0 q|\mathbf{S}\mathbf{v}|, \quad |\boldsymbol{\tau}| \le \alpha_0 q, \quad \boldsymbol{\tau} = \boldsymbol{\tau}^{\mathrm{T}} \quad \text{and} \quad \operatorname{tr}\boldsymbol{\tau} = 0, \quad (56) $$

$$ \operatorname{div}\mathbf{v} = 2\alpha_0|\mathbf{S}\mathbf{v}| - \gamma_0\sqrt{q}, \quad (57) $$

$$ \partial_t\psi + \phi_0\operatorname{div}\mathbf{v} = 0, \quad (58) $$
where $\lambda_0 = \frac{\varepsilon}{\mathfrak{F}\tau^2} = \mathcal{O}(1)$ and $\gamma_0 = \frac{2\alpha_0(\phi_{\max} - \phi_0)}{\mathfrak{D}\mathrm{i}(\phi_0 - \min)} = \mathcal{O}(1)$.

In this system, the evolution of the volume fraction $\psi$ is decoupled from the velocity-stress system (55)–(57). This independent system on $(v, q, \tau)$ is similar to the system (9)–(11) (in which viscosity has been added). In other words, the system studied in the first sections can be seen as an approximation of the complete system with the $\mu(I)$-rheology.

## 4.3 Other models and other possible works

**Case of a granular material immersed in a liquid** Note that in a recent paper by Barker et al [4], the authors use the $\mu(J)$-rheology, instead of $\mu(I)$, to describe fluidised granular flow, i.e. granular material immersed in water. The dimensionless number $J$, defined by $J = \eta_f |\mathbf{S}|/p$ is used for granular flows with a low Stokes number ($St = \rho_0 d^2 |\mathbf{S}|/\eta_f$) and is therefore well suited for granular flows in liquids. As in the case of the models using the $\mu(I)$-rheology discussed in the introduction, instability problems are frequent, and authors often add regularization terms without giving rigorous arguments for the stability of the resulting system. For example, the models introduced in [22, 23] are relatively close to what is proposed here, but the adapted closure (by adding a time derivative on pressure terms) does not seem to ensure stability or the existence of a solution. The study proposed in the present article might be adaptable to such configuration.

**Shallow water model** In gravity flow applications, the domain geometry can be highly anisotropic. Typically, pyroclastic flow lengths are of the order of several kilometers, while their heights do not exceed a few meters. In this context, many authors use models averaged over the transverse component of the flow, thus reducing the number of unknowns and variables in the problem, see for instance [8, 9, 10]. The natural question, then, is what the model introduced in this article might yield when averaged vertically.

**Numerical simulations** The aim of this type of model is to go as far as possible and run numerical simulations to compare with experience. Actually work is in progress to develop a numerical method, in line with the theoretical results presented here, i.e. retaining as many of the proven properties as possible (bounds on volume fraction, dissipated energy, etc.).

## Acknowledgments

This is contribution no. xxx of the ClerVolc program of the International Research Center for Disaster Sciences and Sustainable Development of the University of Clermont Auvergne.

This project was supported by the French institute of Mathematics for Planet Earth (iMPT) and by the project ComplexFlows of the PEPR Math-VivEs, ANR-23-EXMA-0004.

Declaration of Interests: The authors report no conflict of interest.

## References

[1] A. Abbatiello, M. Bulívcek, T. Los, J. Málek, and O. Souvcek. On unsteady flows of pore pressure-activated granular materials. *Zeitschrift für angewandte Mathematik und Physik*, 72(1):6, 2021.

[2] A. Abbatiello, T. Los, J. Málek, and O. Souvcek. Three-dimensional flows of pore pressure-activated bingham fluids. *Mathematical Models and Methods in Applied Sciences*, 29(11):2089–2125, 2019.

[3] B. Andreotti, Y. Forterre, and O. Pouliquen. *Les milieux granulaires-entre fluide et solide: Entre fluide et solide*. EDP sciences, 2012.

[4] T. Barker, J. Gray, D. Schaeffer, and M. Shearer. Well-posedness and ill-posedness of single-phase models for suspensions. *J. Fluid Mech.*, 954:A17, 2023.

[5] T. Barker, D. Schaeffer, M. Shearer, and J. Gray. Well-posed continuum equations for granular flow with compressibility and $\mu(I)$-rheology. *Proc. R. Soc. Lond. A*, 473(2201):20160846, 2017.
[6] T. Barker, D. G. Schaeffer, P. Bohórquez, and J. Gray. Well-posed and ill-posed behaviour of the $\mu(I)$-rheology for granular flow. *J. Fluid Mech.*, 779:794–818, 2015.

[7] A. Beck. *First-order methods in optimization*. SIAM, 2017.

[8] F. Bouchut, E. D. Fernandez-Nieto, A. Mangeney, and G. Narbona-Reina. A two-phase shallow debris flow model with energy balance. *ESAIM: Mathematical Modelling and Numerical Analysis*, 49(1):101–140, 2015.

[9] F. Bouchut, E. D. Fernández-Nieto, A. Mangeney, and G. Narbona-Reina. A two-phase two-layer model for fluidized granular flows with dilatancy effects. *J. Fluid Mech.*, 801:166–221, 2016.

[10] F. Bouchut, E. D. Fernández-Nieto, A. Mangeney, G. Narbona-Reina, et al. Dilatancy in dry granular flows with a compressible $\mu(I)$-rheology. *Journal of Computational Physics*, 429:110013, 2021.

[11] F. Boyer and P. Fabrie. *Mathematical Tools for the Study of the Incompressible Navier-Stokes Equations and Related Models*, volume 183. Springer Science & Business Media, 2012.

[12] E. C. Breard, L. Fullard, J. Dufek, M. Tennenbaum, A. Fernandez Nieves, and J. F. Dietiker. Investigating the rheology of fluidized and non-fluidized gas-particle beds: implications for the dynamics of geophysical flows and substrate entrainment. *Granul. Matter*, 24(1):34, 2022.

[13] H. Brezis. Analyse fonctionnelle. *Théorie et applications*, 1983.

[14] L. Chupin and T. Dubois. Non-isochoric stable granular models taking into account fluidisation by pore gas pressure. *Journal of Fluid Mechanics*, 979:A14, 2024.

[15] L. Chupin, T. Dubois, M. Phan, and O. Roche. Pressure-dependent threshold in a granular flow: Numerical modeling and experimental validation. *J. Non-Newton. Fluid Mech.*, 291:104529, 2021.

[16] L. Chupin and J. Mathe. Existence theorem for homogeneous incompressible navier–stokes equation with variable rheology. *European Journal of Mechanics-B/Fluids*, 61:135–143, 2017.

[17] L. Fang and Z. Guo. Global weak solutions to a three-dimensional compressible non-newtonian fluid. *Communications in Mathematical Sciences*, 20(6):1703–1733, 2022.

[18] K. Hutter and K. Rajagopal. On flows of granular materials. *Continuum Mechanics and Thermodynamics*, 6:81–139, 1994.

[19] P. Jop, Y. Forterre, and O. Pouliquen. A constitutive law for dense granular flows. *Nature*, 441:727–730, 2006.

[20] J. Malek, J. Necas, M. Rokyta, and M. Ruzicka. *Weak and Measure-Valued Solutions to Evolutionary PDEs*. Chapman and Hall, 1996.

[21] G. MiDi. On dense granular flows. *Eur. Phys. J. E*, 14:341–365, 2004.

[22] E. Montellà, J. Chauchat, C. Bonamy, D. Weij, G. Keetels, and T. Hsu. Numerical investigation of mode failures in submerged granular columns, flow, 3, e28, 2023.

[23] E. P. Montellà, J. Chauchat, B. Chareyre, C. Bonamy, and T.-J. Hsu. A two-fluid model for immersed granular avalanches with dilatancy effects. *Journal of Fluid Mechanics*, 925:A13, 2021.

[24] F. Radjai and F. Dubois. *Discrete-element modeling of granular materials*. Wiley-Iste, 2011.

[25] J. A. Robinson, D. J. Holland, and L. Fullard. Complex behavior in compressible nonisochoric granular flows. *Phys. Rev. Fluid*, 8(1):014304, 2023.

[26] O. Roche. Depositional processes and gas pore pressure in pyroclastic flows: an experimental perspective. *Bull. Volcanol.*, 74(8):1807–1820, 2012.

[27] S. Roux and F. Radjaï. Texture-dependent rigid-plastic behavior. In *Physics of dry granular media*, pages 229–236. Springer, 1998.
[28] D. Schaeffer, T. Barker, D. Tsuji, P. Gremaud, M. Shearer, and J. Gray. Constitutive relations for compressible granular flow in the inertial regime. *J. Fluid Mech.*, 874:926–951, 2019.

[29] V. Zhikov. On the weak convergence of fluxes to a flux. In *Doklady Mathematics*, volume 81-1, 2010.