{'fraction_non_alphanumeric': 0.04799183117767188, 'fraction_numerical': 0.014068527342863625, 'mean_word_length': 4.1441494018091625, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
FLUIDISATION BY PORE PRESSURE OF DENSE GRANULAR FLOWS: NUMERICAL SIMULATIONS VERSUS EXPERIMENTS

. A model accounting for ﬂuidisation by pore gas pressure in dense granular ﬂows is presented. A viscoplastic rheology, based on the Drucker-Prager criterium, is used to describe the granular medium which is a mixture of air and glass beads. The pore gas pressure, which satisﬁes an advection-diﬀusion equation, reduces the friction between the particles and thus the value of the apparent viscosity. As a consequence, dense ﬂuidised granular ﬂows can travel longer distances. In laboratory experiments, the run-out distance reached by dense granular columns when collapsing is almost doubled when ﬂuidisation is applied. This fundamental result, in the context of pyroclastic density currents, is reproduced by numerical simulations performed with the ﬂuidised model.


INTRODUCTION
Pyroclastic density currents (PDCs) are one of the most important hazards occurring in volcanic eruptions. Indeed, PDCs are mixtures of solid particles (pyroclasts and lithic fragments) and air that are able to travel over long distances from the eruption source, which may exceed 100 km in some cases. PDCs are generally induced by the collapse of a volcanic dome or an eruptive column. They often consist of a dense basal flow overlain by a dilute turbulent flow. The dense part behaves like a fluid and is able to travel long distances. Understanding the mechanisms responsible for this peculiar behaviour of concentrated PDCs is one of the major scientific question related to volcanic processes. A possible answer is a fluidisation effect of the pore gas pressure that reduces friction between particles. This question has been addressed through laboratory experiments in [1]. The main result is that the run-out distance reached by a dense granular columns when released is almost doubled when fluidisation is applied. Numerical simulations with an averaged (i.e. one layer) model has been used in [2] to reproduce experiments.
In the case of the collapse of eruptive columns, the initial height of the collapsing granular mixture is significantly greater than the width so that averaged models could not, in principle, be applied. In this context, we introduce a non-averaged model, which also accounts for fluidisation by pore gas pressure. The model is based on the mass and momentum conservation equations with a viscoplastic rheology. The yield stress is a Drucker-Prager criterium depending on the solid (inter-particle) pressure. As in [3], the pore gas pressure obeys a diffusion-advection equation, which is obtained by assuming that a Darcy law is satisfied in the solid-fluid mixture and that the gas is perfect. In our work, the air surrounding the granular deposit during its slump is also taken into account by considering the incompressible Navier-Stokes equations. In order to obtain a unified model and unknowns (velocity, pressure) for both the granular and fluid phases a level set formulation is used as in [4].
The viscoplastic rheologies have the particularity of being defined with a non-differentiable term which induces mathematical and numerical difficulties. As in [5,6], we overcome these difficulties by rewriting the definition of the plastic part of the stress tensor as a local projection. Based on this formulation, a bi-projection scheme for the time discretisation of the Bingham equations has been proposed and analysed in [5] and extended in [6] to viscoplastic fluids with variable density, viscosity and yield stress. In the bi-projection scheme, the plastic part of the stress is computed with a Picard fixed point procedure accelerated by adding a pseudo-relaxation term which makes the convergence geometric. This approach avoids the use of regularisation methods [7] and the formalism of variational inequalities [8]. The performance of the bi-projection scheme was evaluated by simulating Bingham flows in a lid-driven cavity in [5] as well as by successfully reproducing experiments of the collapse of non fluidised dense granular columns, i.e. without pore gas pressure, in [4]. The model, based on the µ(I)-rheology [9], reproduces the internal dynamics of such granular flows characterised by the presence of an interface which separates a static basal deposit and an upper moving part. The interface migrates upwards in the course of time and controls the deposition rate of the particles.
The aim of this paper is to reinforce the results obtained in [10] on relatively coarse grids, which show that incorporating the effect of pore gas pressure in the model recovers a fundamental result observed experimentally: fluidisation allows granular flows to propagate over longer distances from their source compared to non fluidised flows. A grid twice the size in [10] is used in the work presented here. Moreover, instead of using a space and time variable granular viscosity as in [10], we have chosen to take a constant value as in a previous work on nonfluidised columns [4]. In [10], emphasis was given on the analysis of simulation results in many different configurations while here we focus on the description of the mathematical model and its numerical approximation.
The paper is organised as follows. In Section 2, the mathematical model is described. In Section 3, the numerical schemes are summarised. In Section 4, the results of a numerical simulation of the collapse of a fluidised column of glass beads, with an aspect ratio of 2 over a horizontal plane, are presented and compared with experimental results. The granular phase (blue), which is a solid-fluid mixture, occupies the region Ω g (t) while in the remaining part Ω f (t) = Ω\Ω g (t) only air is present. The interface ∂Ω g (t)∩∂Ω f (t) between the granular medium and the air is denoted Γ gf (t).

The governing equations
As the experimental apparatus used in [1] is a narrow channel (10 cm wide and 3 m long), we assume that the granular flow is mainly two-dimensional. In that context, the computational domain is Ω = (0, L)×(0, H). At time t ≥ 0, the granular medium, consisting of 58% per volume unit of glass beads surrounded by air, occupies the subdomain Ω g (t) while Ω f (t) = Ω\Ω g (t) contains only air (see Figure 1). Therefore, the volume concentration of particles φ is equal to 0.58 in Ω g (t) and 0 in Ω f (t).
In Ω f (t), we assume that the flow is governed by the incompressible Navier-Stokes equations, namely the velocity field u f and the pressure p f are solutions of where η f = 2 × 10 −5 Pa s is the air viscosity, ρ f = 1 kg m −3 is the air density and g = (0, −g) with g the gravitational constant.
In Ω g (t), the particles are spherical glass beads of density ρ s = 2500 kg m −3 and mean size d s = 75 µm. As this study is restricted to dense granular flow, the volume concentration φ is assumed to be constant during the slump of the granular mass, namely φ = 0.58, so that the velocity field u s of the granular medium is solenoidal. The granular medium is modelled as an incompressible fluid with a specific rheology. Due to the presence of air between the glass beads (see Figure 1), the total granular pressure is decomposed as the sum of the solid (effective) pressure p s and the pore gas pressure p f . The dynamics of the latter is modelled by a diffusionadvection equation. This approach is supported by experimental results (see [11,1]). Hence, the velocity field u s , the pressures p s and p f are solutions of the following set of equations where κ f is a diffusion coefficient. Note that, as in [3], the pore gas pressure equation is obtained from the mass conservation equation of the fluid phase in Ω g (t) by assuming that the density ρ f is proportional to the pore pressure (perfect gas law) and that the fluid velocity is given by the Darcy law.

The rheology of the granular medium
A granular medium flows only if the stress exceeds a threshold value, the yield stress. Otherwise, it does not deform and behaves like a solid. To date, the most advanced model for this particular behaviour is the µ(I)-rheology [9] which defines the deviatoric stress tensor T s by where α is the internal friction angle of the particles, µ ∞ is an asymptotic friction coefficient and I 0 is a dimensionless number. The inertial number I is defined by In (3) and (4), D(u s ) = 1 2 (∇u s + T ∇u s ) is the strain-rate tensor. Also, the tensorial norm |λ| = 1 2 i,j λ 2 ij , for any tensor λ ∈ R 2×2 , is used. By reporting the expression (4) in (3), we easily deduce the following form of the deviatoric stress in the zones where deformation occurs (i.e. D(u s ) = 0), with η s (|D(u s )|, p s ) = (µ∞−tan(α))ps 2|D(us)|+ I 0 ds ps ρs . The above expression of the deviatoric stress tensor modelled by the µ(I)−rheology separates the diffusion and plastic terms; the latter satisfy a Drucker-Prager yield criterium. Note that by replacing the time and space varying viscosity η s (|D(u s )|, p s ) in (5) by a constant value η s , we obtain the Drucker-Prager viscous rheology, hereafter referred to as the D-P model. As in [5], we introduce the symmetric traceless tensor Σ defined by so that (5) rewrites The tensor Σ can be equivalently defined as a projection. Indeed, we have where P is the orthogonal projection on the unit ball in the space of symmetric traceless tensor (see [5] or [6] for the details). Note that this formalism avoids the use of a regularised viscosity (see [7] for instance). Indeed, the rheology for a yield-stress fluid leads to an infinite viscosity in regions where the strain-rate is zero. Regularisation methods replace the singular viscosity by a large but finite one. This is not necessary with the formalism based on a projection introduced above: the effective viscosity is now everywhere well-defined. The model for the granular flow is now complete and reads φρ s ∂ t u s + u s · ∇u s + ∇(p s + p f ) = φρ s g + div(2η s D(u s )) + div(tan(α) p s Σ), div u s = 0, (7b)

Boundary conditions
The equations (1) and (7) must be supplemented with appropriate boundary conditions. Free, Navier and Dirichlet boundary conditions are imposed to the air velocity field u f on the boundary ∂Ω ∩ ∂Ω f (t) as shown on Figure 2 (in the particular case where t = 0). As in [12], Coulomb friction conditions are imposed to the granular flow on the boundary ∂Ω ∩ ∂Ω g (t), namely they read where α b is the basal friction angle, σ s,t (resp. σ s,n ) is the tangential (resp. normal) component of the total stress tensor σ s = T s − 1 2 p s I and u s,t is the tangential component of the granular velocity field. By analogy with (6), we introduce a vector s with Euclidean norm less or equal than unity, defined by s = − us,t |us,t| , if u s,t = 0, |s| ≤ 1, if u s,t = 0, so that (8) rewrites : σ s,t = tan(α b )(−σ s,n ) + s. As previously, s can be equivalently rewritten as a projection on the unit ball in R 2 , namely During the experiment a constant air flux is applied at the bottom of the granular column in the reservoir, that is for x ∈ (0, L c ) × {0} (see Figure 2). This is modelled by applying in the numerical simulations a constant pressure gradient ∂p f ∂n = −0.9φρ s g on this part of the boundary. The factor 0.9 is based on experimental constraints [13]. Everywhere else on ∂Ω homogeneous Neumann condition is applied. Accordingly, the pore gas pressure at t = 0 is defined by p f (x, y, t = 0) = 0.9φρ s g(H c − y) in Ω g (t = 0) = (0, L c ) × (0, H c ) and 0 in Ω f (t = 0).

Level set formulation
The equations (1) and (7) describing the evolution of the flow in Ω f (t) and Ω g (t) can be reformulated in one set of equations on the whole computational domain Ω by using a level-set formulation, namely The level-set function Φ characterises the interface between the air and the granular medium. Indeed, at each time t > 0 , we define: Γ gf (t) = {x ∈ Ω; Φ(x, t) = 0}, Ω g (t) = {x ∈ Ω; Φ(x, t) > 0} and Ω f (t) = {x ∈ Ω; Φ(x, t) < 0}. At time t = 0, the level-set function is initialised as the signed distance to the interface. In order to keep Φ as close as possible to the signed distance function, it must be periodically reinitialised in time; this is achieved by applying the redistancing algorithm proposed in [14]. In (10b), the viscosity η(Φ) is defined as the harmonic mean of η s and η f (see [4] for details). Let a parameter > 0, we introduce the following regularised Heaviside function The density and the diffusion coefficient are defined with the help of H by The parameter κ ∞ is set to a large value (10 16 in practice), so that the solution of the pore gas pressure equation (10f) is a smooth extension by zero outside Ω g (t) of the solution of (7d). The jumps across the interface Γ gf (t) of the yield stress and the friction coefficient are also regularised by using the function H , indeed 3 The numerical schemes

A bi-projection time scheme
Let δt > 0 a time step and t n = nδt a sequence of associated discrete times. Let us assume that, at t n , the approximations (u n , p n , p n f , Φ n , Σ n , s n ) of the solutions (u, p, p f , Φ, Σ, s) of the set of equations (10a)-(10f) are known and that u n is solenoidal. The time discretised scheme consists in the steps described below.
The redistancing algorithm [14] is applied about every ten time steps.

The spatial discretisation
We denote by (u, v) the components of the velocity field u. The computational domain Ω = (0, L)×(0, H) is discretised by using a Cartesian uniform mesh. This results in a sequence of nodes denoted (x i , y j ) for a given mesh size. Let a computational cell K ij = (x i , x i+1 )×(y j , y j+1 ). The discrete unknowns of the velocity (resp. pressure) are placed in the middle of the edges (resp. in the centre) of the K ij mesh. The discrete unknowns associated with Σ and p f are also placed at the centre of the mesh cell. Figure 3 summarises this staggered arrangement of the unknowns. Cell-centred second-order finite volume schemes are applied to discretise first and second-order partial derivatives with respect to the spatial directions in equations (10b), (10c) and (10f). The transport equation of the level set function (10a) is discretised with an upwind WENO5 scheme. The code is parallel, written in F90 and using the MPI library to handle communications between processes. The PETSc library is used to assemble and solve the linear systems.

Numerical simulations versus experiment
We consider the collapse of a dense column of glass beads of width L c = 20 cm and height H c = 40 cm initially at rest and fluidised in a reservoir which is connected with a horizontal channel with a smooth base through a gate opened at t = 0 s. In the simulation, we assume that the whole granular column is released at t = 0 s. The internal friction angle of the material is α = 27°and the basal friction angle at the bottom and on the vertical back-wall of the channel is α b = 15°(see [4]). In agreement with the experiment, the pore pressure at time t = 0 is set to 90% of the weight of the particles, that is : p f (x, y, t = 0) = 0.9 φρ s (H c − y) in Ω g (t = 0) and p f (x, y, t = 0) = 0 in Ω f (t = 0). Results of a simulation with a mesh size h = 1.7578 mm, corresponding to 256 cells in the vertical direction, and with κ f = 0.035 m 2 s −1 are reported ( Figure 4). As in [4], the D-P model is used as rheology for the granular flow with an apparent viscosity η s = 0.1 Pa s. A simulation with the same initial configuration but on a two times coarser mesh and with the variable granular viscosity of the µ(I)-rheology has been reported in [10]. The value of the time step, δt = 10 −5 s, is chosen so that the CFL number remains A. Aravena, L. Chupin, T. Dubois and O. Roche The model underestimates by about 8.5% the run-out distance found in the experiment (see Figure 5). The simulation predicts a distance of 1.6874 m from the reservoir gate, whereas 1.8443 m was found in the experiment. The run-out distance is almost doubled when fluidisation is applied (see Figure 5). Note that in [10] the granular flow profile extended only to 1.4 m from the reservoir gate which corresponds to an underestimate of about 32%. The simulation (see Figure 5) reproduces reasonably well the three phases of propagation described in [1] and the duration of each phase. Note that the granular flow stops a bit earlier in the experiment. As shown on Figure 4, the height of the deposit at the back wall of the channel rapidly decreases (t ≤ 0.7 s) and then becomes stationary. During the initial phase, the height of the deposit is underestimated whereas a perfect match is obtained later. The thickness of the deposit is slightly underestimated for x − L c ∈ [0.1, 0.5] m and overestimated at longer distance from the gate. However, the overall profile of the granular flow deposit is well predicted with a peak of the thickness near the reservoir x − L c = 0 m and a gentle decrease downstream. By comparing Figure 5 at t = 1.2 s and Figure 2f in [10], we observe that a better approximation of the profile is obtained with the D-P model. Especially, near the front position, namely for x − L c ≥ 1.2 m the profile is thicker when the granular viscosity is kept constant. This remark should be mitigated by the fact that in the present work the grid is two times as small. A fair comparison of the models should be done with the same numerical parameters.
As a summary, in terms of run-out distance, dynamics of the collapse and shape of the granular mass, the model provides results in good agreement with the experimental measurements. The effect of the pore gas pressure observed in experiments is reproduced by the model.