An Extension of the Drift-Flux Model for Submarine Granular Flows

The flow of granular material submerged in water is a regularly occurring phenomena. This is the case, for example, in breaching [1, 2], high speed erosion [3], underwater avalanches [4, 5] and debris flows [6].

In these phenomena we have to deal with a large range of different particle concentrations, ranging from pure water (0%) up to dense flows where the friction between particles dominates (~55%). In the friction dominated regime, the mixture does not act as a Newtonian fluid, and the standard Navier–Stokes equations are therefore no longer valid. Instead the material remains rigid, unless the shear stresses inside the material are high enough. The necessary shear stress for motion is dependant on the pressure between sand particles, also known as the dispersed phase pressure. The higher this pressure, the larger the shear stress needs to be before movement is initiated.

Many authors simplify this behaviour to that of a fluid with a yield stress. One of the first known attempts to represent dense granular material as a fluid with a shear stress is by Lalli and co-workers [7, 8]. They use the drift-flux model [9], meant for flows with a concentration up to 30%, as their starting point. They extend the model by introducing a yield stress for higher concentrations, where friction between particles introduces a yield stress. To calculate the dispersed phase pressure, necessary to calculate the yield stress, they assume that the pressure is hydrostatic. This assumption is not valid near steep slopes. For the flow after yielding they assume a constant ratio between the shear rate and the dispersed phase pressure. This is only a first assumption and not based on any experiments.

A French research group, GDR MiDi [10], investigated the results of experiments and discrete particle simulations of different types of dense granular flows in air. Based on these results they proposed a relation between the normal pressure and the shear rate, after yielding of the material. This relation for shear stress was implemented in a 3D numerical simulation by Jop et al. [11]. This model is similar to that of Lalli and co-workers [7, 8], but includes an experimentally derived relation between normal pressure and shear rate after yielding, instead of the fixed ratio proposed by Lalli and co-workers.

The framework was extended to granular flows submerged in water by Cassar et al. [12], using the work of Courrech du Pont [13]. They assume a hydrostatic pressure to calculate the dispersed phase pressure. Pailha & Pouliquen [5] extended the model by including the effects of dilatancy and pore-pressure feedback, in a vertical 1D model, to investigate the initiation of avalanches. They also assume hydrostatic pressure to calculate the dispersed phase pressure.

Bohorquez [14] created a three-phase mixture model, to model the flow of dense suspensions in air. His model did not calculate the dispersed phase pressure. Iverson [6] also created a model which includes all the aforementioned phenomena, to model debris flows. He also included the elasticity of the grain skeleton [15] attempted to model the collapse of rectangular granular piles in water. In the process the dilatancy effects play an important role. In his model, [15] modelled the effects of dilatancy by temporarily increasing the viscosity of the sand. This approach is limited to the collapse of rectangular piles. It also does not take the physics behind the process into account.

All the aforementioned models either calculate the dispersed phase pressure directly from the balance equations, or assume a hydrostatic pressure distribution. In the case of our multiphase approach to modelling, it is not possible to calculate this pressure directly from the balance equations. It is possible to assume a hydrostatic balance, but this assumption looses validity near steep slopes. Therefore, we propose a new way to calculate the dispersed phase pressure.

In our study we choose to use the drift-flux model as our starting point. The drift-flux model is derived starting from the balance equations for multiphase flow. Here we use to multiphase formulae derived by [16].

We simplify these equations by summing the continuity and momentum equations, giving us mixture continuity and mixture momentum equations. Here the drift-flux model developed by Manninen et al. [17] is explained.

The continuity equation is summed over all phases. We obtain

where $\rho_m$ is the mixture density, and $u_m$ is the velocity of the centre of gravity of the mixture. These are defined as

where $\underline{\boldsymbol{u}}_k$ is the velocity of phase k , and $\rho_k$ the density of phase k. The momentum equation is summed as well. This leads to the following equation

The term $\underline{\boldsymbol{M}}_{\mathrm{m}}$ is the influence of the surface tension force on the mixture, which can be taken as zero in the case of solid particles. The term $p_m$ is the mixture pressure, and is equal to the linear combination of the dispersed phase and continuous phase pressures. The term $\underline{\underline{\tau}}_m^{\mathrm{D}}$ represents the diffusion stress due to the differences between phase velocities and mixture velocity.

The phase velocities, $\underline{u}_k$, are not calculated in this model. Therefore a closure relation is needed. Usually this is accomplished by rewriting $\underline{\underline{\tau}}_m^{\mathrm{D}}$ in terms of the drift-flux of the dispersed phase, $v_{d j}$. In the drift-flux model it is assumed that a local equilibrium drift-flux will be established over relatively short length scales. This assumption is correct for small particles $(d<200 \mu \mathrm{~m})$ [17]. The drift-flux velocity is the difference between the velocity of a phase and the volumetric flux of the mixture. The volumetric flux, $\boldsymbol{j}$, is the velocity of the volume centre and is defined as

where $\alpha_k$ is the volume concentration of phase $k$.

The drift flux of phase $k$ is defined as

To close the system we compute the drift-flux of the dispersed phase with the formula of Richardson and Zaki [18]

Finally the continuity equations of the dispersed phase, $d$, is rewritten using the mixture velocity, $\underline{u}_{m^{\prime}}$, and the drift-flux velocity, $\underline{v}_{d j^{\prime}}$

The mixture model is sufficient for modelling flows with concentrations up to 30% [19], but for the bed the effects of friction between particles have to be taken into account. For this we propose to model the sand bed as a Bingham fluid as first proposed by Lalli and Di Mascio [7]. In the Bingham rheology the shear stress tensor, $\underset{=}{\tau}$, is defined as.

Here $\mu_c$ is the dynamic viscosity of water in $\mathrm{Pa} \mathrm{s} \dot{\gamma}$ is the second invariant of the deformation tensor, $\underline{\underline{\gamma}}$, and, using the Einstein summation convention, is defined as follows:

and

When the shear stress is below a given yield stress, $\tau_y$, the rate of strain, characterised by $\dot{\gamma}$, is 0 , meaning no flow. When the shear stress is larger than the yield stress, the stress tensor is given by Eq. (9). Lalli and Di Mascio [7] applied the yield stress of Coulomb:

where $p_{\mathrm{d}}$ is the dispersed phase pressure, $\phi$ is the internal friction angle, and $c$ is the cohesion of the soil.

The Bingham framework cannot be applied directly in a numerical method. Usually the Bingham fluid is applied by adjusting the viscosity:

To avoid division by zero problems, a small regularization parameter is added. Other approaches were reviewed by [20].

GDR MiDi [10] investigated the results of experiments and discrete particle simulations of different types of dense granular flows in air. Based on these results they introduced the inertial number, $I$. I is the ratio between the macro time scale of granular flow $(|\dot{\gamma}|)$, and the micro time scale $\left(\sqrt{p / \rho_d / d}\right)$. Thus, $I$ is defined as

Because $I$ is a local parameter it cannot take into account non-local effects such as arching or grain clusters. However, this rheology seems to be able to accurately predict many different types of granular flows [11, 12, 21-23].

Jop et al. [11] created a model for flow of granular material based on the results of GDR MiDi [10] for dry granular flow. They calculate the yield stress as follows:

where the friction coefficient $\mu(I)$ is:

Cassar et al. [12] investigated submarine flows of granular materials down a rough incline. In this case the interstitial fluid plays an important role in the flow. They found that the flow can still be described by the μ(I)-rheology if a micro time scale is chosen which takes the effect of this interstitial fluid into account. The alternative time scale is based on the work of Courrech du Pont et al. [13]. They replaced the inertial number with,

For higher Stokes numbers $\left(S t=\left(1 / 18 \sqrt{2} \sqrt{\left.\rho_d\left(\rho_d-\rho_c\right) g \sin \phi\right)} / \mu_c\right)\right.$, inertia effects become important. The inertial number in this regime is.

Where $C_d$ is the draf coefficient of a sand particle.

The dispersed phase pressure term, $p_{\mathrm{d}}$, in Eq. (15) follows directly from the balance equations for dry granular flows [11]. However, for submerged granular flows it does not.

A common way to calculate $p_d$ is to assume hydrostatic pressure [5, 12, 15]. In this case the velocities are assumed 0, or constant in time and space. This reduces the momentum equations to

Without neglecting any terms, the gradient of the dispersed phase pressure is

The difference between the actual and estimated pressure gradient, in the case of the hydro-static pressure assumptions is thus

This is a reasonable assumption when working with flat or almost flat beds, with small accel-erations. But we want to do calculations with non-stationary flows, with large accelerations. Even in slowly flowing granular material a large error is caused by the viscous shear stress term of the dispersed phase.

The effective viscosity, $\mu_{\mathrm{e}} \mathrm{ff}$, can go up to very large values, especially in quasi-static regions. This causes the shear stress of the dispersed phase to have significant values, even when the velocity gradients seem negligible. Therefore we propose to calculate the dispersed phase pressure in a different way, that includes all terms of Eq. (21).

The mixture pressure is already computed by the model, and includes the effects of the dispersed phase shear stress. We can use this fact to calculate a better approximation of the dispersed phase pressure gradient

As a first approximation we can assume hydrostatic water pressure, giving us

where $p_d^*$ is the estimated dispersed phase pressure. The difference between to approximation and the real gradient is now

We can plug this $\boldsymbol{u}_{\boldsymbol{c}}$ into the continuous phase momentum equation and calculate $\nabla p c$. We then use Eq. (24) to calculate $\nabla p d$. With this method the calculated $\nabla p_d$ will be equal to that of Eq. (21). It is important to note that in all cases the dispersed phase pressure is assumed to be isotropic.

We extended the mixture formula with the effective Bingham viscosities using Eqs. (12), (13) and (17). To test the model we started with heaps of sand with initial angles of 45 and 30 degrees, and an angle of repose of 35 . The sand was given a density of $2,650 \mathrm{~kg} \mathrm{~m}^{-3}$ and an initial porosity of 0.45 . The results are shown in Figure 1.

Figure 1a shows results of the simulation with a starting angle of 45 degrees. The contours of the heap are shown at the start of the simulation, and with intervals of 1 second until 5 seconds. After 5 s the sand has come to rest it can be seen that the sand heap came at rest at an angle slightly below the internal angle of friction of 35 .

Figure 1b shows results of a similar simulation but with an initial angle of 30 degrees. It can be seen that in this case there is only some slight erosion at the toes of the heap and at the top.

The results show that we are able to recover the internal fraction angle with the model. When the initial angle of the slope is higher than the internal angle of friction, the soil is unstable and starts to flow. The slope starts to reduce and when the angle is below the internal friction angle the flow of material is slowed down. This leads to a heap with an angle just below the internal angle of friction.

When the initial angle of the heap is below the internal friction angle, the heap is stable as is expected. However, we do see a little flow of material at certain points. This can be due to our assumption that the dispersed phase pressure is isotropic, which is not the case. Also our use of a Cartesian grid could introduce errors, since this leads to a ‘staircase’ edge of the heap, instead of a straight boundary.

The time it takes for the collapse of the heap to complete, around 3 seconds, is in the same order of magnitude as experiments in which collapse of a granular column in water, without pore pressure effects, is investigated [24–26]. This indicates that the effective viscosity and dynamics calculated by the model is also in the right order of magnitude.

We have created a model that is able to model both fluid as well as dense granular materials. Compared to similar models we have devised a new method of calculating the dispersed phase pressure. This new method allows us to more accurately determine the dispersed phase pressure. This is particularly useful in cases with steep slopes, where the use of a hydrostatic pressure seems to fail. With this model we are able to reproduce the stability of sand piles at their internal angle of friction, with time scales in the right order of magnitude.

There are, however, also certain aspects of dense granular materials that the model is not able to capture. The dispersed phase pressure is assumed to be isotropic, which is not the case in reality. But we believe that in our case the anisotropy is limited, and thus has a limited influence on the results. We can also not capture non-local effects such as arching or grain clusters.

However, we believe that our model is able to model is able to correctly capture the behaviour of dense granular material. It is not advised to use the model for detailed soil mechanics, such as soil stability analysis, but it provides a good alternative when less detail is required.