# IJEGE-11_BS-Minatti-&-Pasculli

*Italian Journal of Engineering Geology and Environment - Book www.ijege.uniroma1.it © 2011 Casa Editrice Università La Sapienza*

*DOI: 10.4408/IJEGE.2011-03.B-052*

**SPH NUMERICAL APPROACH IN MODELLING 2D MUDDY DEBRIS FLOW**

*alii*, 2001). The most widely used approaches to model

debris flows are usually related either to “continuum”

or to “granular” mechanics. A combination of them can

sometimes be used, depending on the type of debris

flow under study. The content of this paper belongs to

the former kind of approach where the Navier-Stokes

equations are considered as the master tool. In the

“continuum” framework, a quite general mathematical

model, that can be used to model a physical phenom-

enon like the one discussed in this paper, would consist

of two sets of equations: a first set for the liquid phase

and a second one for the solid phase. Interactions terms

can then be added to the equations in order to model

the forces exchanged by the two phases. A model of

this kind can be derived from the mixture theory (a

*et alii*, 1976): an implementation of such theory

logical models for each phase is suggested.

in the study of water-solid mixtures have shown that

the clay fraction can be of great importance too, as

it influences the granular interactions. In particular, it

is reported by l

*et alii*(1997) that debris flow

their two-phase nature, behave like viscoplastic fluids

and that a single phase model of such kind can repro-

duce their physics to a good degree of approximation.

We then restrict our analysis to the specific case of

**ABSTRACT**

sets of equations related to the fluid and solid phases,

as considered by the debris flow mixture theory, have

been simplified in only one set of equations, consider-

ing just one equivalent material. Then the Herschel-

Bulkley fluid constitutive equations have been select-

ed. The correct parameters of the H

the behaviour of mudflows. The final mathematical

model, has been solved numerically with the smoothed

particle hydrodynamics (SPH) method. SPH is a par-

ticle mesh-free Lagrangian method, well suitable for

computing highly transitory free surface flows of com-

plex fluids in complex geometries. Finally a laboratory

experimental test has been selected for comparison.

Satisfactory results have been achieved. Nevertheless,

further parametric analyses will be carried out and fur-

ther considerations about both constitutive equations

and numerical improvements will be employed and

discussed in future papers.

**K**

**ey**

**w**

**ords**

**:**SPH, 2D numerical modelling, muddy debris flow**INTRODUCTION**

the effect of gravity (P

*et alii.*1987; t

*et alii*, 1996; i

*et*

*L. MINATTI & A. PASCULLI*

*5th International Conference on Debris-Flow Hazards Mitigation: Mechanics, Prediction and Assessment Padua, Italy - 14-17 June 2011*

*η*is the local dynamic viscosity of the fluid;

τ

*k*is called the liquid consistency;

*n*is called the power law index;

**D**is the strain rate tensor of eqn. (4).

*k*, and

*n*represent the constitutive law param-

ments (see Tab. 1) in order to correctly reproduce the be-

haviour of a mudflow. The physical meaning of the yield

stress is immediate, representing the stress threshold

below which the fluid starts to behave like a rigid body.

**NUMERICAL MODELLING**

approach. The related bibliography is very extensive and,

for the sake of completeness, we mention the fields of

Computational Fluid Dynamics (CFD) and of Computa-

tional Granular Dynamics (CGD). CFD (C

related to continuum mechanics and granular mechan-

ics. Molecular Dynamics is a further numerical tool

used to study the interactions among particles, exploit-

ing essentially Newton law related to both their transla-

tions and rotations. Nevertheless, according to the scope

of the present paper, an appropriate numerical method,

among the ones suitable for solving the Navier-Stokes

equations, should be selected. In these regards, many ap-

proaches have been proposed: in particular, techniques

capable of handling free surfaces and violent mass fluxes

(mudflows) and simplify the problem by employing a

single equivalent phase model of a viscoplastic fluid.

**GOVERNING EQUATIONS**

motion of an incompressible non-Newtonian fluid. In

order to develop a theoretical framework for a numeri-

cal model we resort to the basic principles governing

the motion of a continuum, namely the mass conserva-

tion and the momentum equations. They can be written

in Lagrangian form as:

*n*

*ρ*is the local density of the continuum

*f*is the body force per unit of mass exerted on the con-

tinuum;

**σ**is the local total stress tensor;

sor decomposition is indicated as follows in the paper:

*p*is the isotropic pressure;

**I**is the unit tensor;

**τ**is the deviatory part of the total stress tensor;

tensor is usually expressed as a function – namely the

constitutive law - of the strain rate tensor

**D**:

Bulkley law, which has also been recently used by

l

*et alii*(2007):

*Tab. 1 - Rheological characteristics and experimental pa-*

*rameters in l*

*AiGle*

*et alii (1994) and l*

*AiGle*

*(1997)*

**SPH NUMERICAL APPROACH IN MODELLING 2D MUDDY DEBRIS FLOW**

*Italian Journal of Engineering Geology and Environment - Book www.ijege.uniroma1.it © 2011 Casa Editrice Università La Sapienza*

ian nature make it very attractive for solving fluid flow

problems where free surface boundary conditions and

large strain rates are involved. The computational do-

main is filled with particles carrying flow field infor-

mation (e.g. pressure, velocity, density) and capable of

moving in space.

a mesh to calculate spatial derivatives is not needed.

can be extended to a 3D case with little effort. The key

idea on which the method is based is the well-known

use of a convolution integral with a Dirac delta func-

tion to reproduce a generic function

*f(x)*:

*w*(it ‘mimics’ the

Dirac delta function), and the generic function

*f(x)*

is reproduced with a convolution integral, which in

a discrete framework takes the form of a summation

over particles:

*x*

*i*

*x*

*j*

*represent the*

*i*and

*j*particle positions in the

*∆A*

*j*

*j*;

Summation is extended to all the particles located

within the support domain of particle

*i.*

*Ω*

*x*

*h*, named

*smoothing length*. The kernel function is zero outside

the support domain and the smoothing length serves

as a scaling parameter for its arguments. It also has

the property of converging to the Dirac function as the

smoothing length approaches to zero.

convolution integral of the function with the kernel

and by using the Gauss-Green formula:

problems.

mon used methods (C

step to be performed. And as far as grid generation is

concerned, the possibilities are numerous: structured

grids, unstructured grids and adaptive meshes.

Navier-Stokes equations involve incompressibility and

convective terms. In order to avoid the latter one, a La-

grangian approach is commonly selected instead of an

Eulerian one.

tions. This is the main reason why a frequent (time

consuming) update of the grid, aimed at lowering the

excessive mesh distortion due to large deformations,

is often necessary.

convective terms, lower the grid generation time and

that are capable of easily taking free surfaces into ac-

count, are very desirable. The “Meshless Finite Ele-

ment Methods” include many of the above necessary

discussed features (see i

*et alii,*2003; 2004;

namics” (SPH), briefly described in the following par-

agraphs and selected in our approach is among these.

A similar tool, the “Particle Finite Element Method”

(PFEM) (o

*et alii*, 2004) is very promising as it

of the difficulties of the SPH method (in particular, the

imposition of boundary conditions), even if it requires

some time, (which can still be reduced by using an

extended Delaunay method), to update grid meshes.

Elements” and the “Alternate Lagrangian Eulerian”

method (ALE) (C

*et alii,*2001).

**OVERVIEW OF THE SPH METHOD**

veloped during the 1970s to solve astrophysical prob-

lems (m

*L. MINATTI & A. PASCULLI*

*5th International Conference on Debris-Flow Hazards Mitigation: Mechanics, Prediction and Assessment Padua, Italy - 14-17 June 2011*

**SPH DISCRETIZATION TECHNIQUE**

*MAIN EQUATIONS*

Bulkley fluid flow governing equations in order to

create an SPH algorithm.

be effectively discretized into SPH equations. Among

the many of them, we indicate m

*et alii*(2003) & t

*et alii*(1994).

*m*represents the particle mass;

*v*

*ij*

*v*

*i*

*v*

*j*

*r*

*ij*

*r*

*i*

*- r*

*j*

*η*

*ij*

*η*

*ij*

*= η*

*ji*

*.*

in the summations is used when calculating the kernel

gradient (e.g.

*hij=(hi+hj)/2*). If this is the case, it can

be shown that particles exchange equal and opposite

forces, making eqn. (12) capable of conserving the

linear momentum of the particles system.

the kernel is differentiated with respect to the

*x’*co-

ordinate;

*n*represents the normal to the support domain bound-

aries (pointing outward).

domain boundaries, as the kernel is zero on the sup-

port domain boundaries. Another case when the term

can be zero is when the support domain is truncated

by the computational domain boundaries but there ex-

ists a boundary condition forcing the function

*f(x)*to

vanish on the boundaries (it may be the case when

*f(x)*

represents a velocity and a no-slip condition has to be

enforced on the computational domain boundaries). If

the first term of the RHS of eqn. (9) is zero, then the

SPH approximation of

*f(x)*gradient takes the follow-

ing form in a discrete framework:

it is possible to offset the errors induced by neglect-

ing the term by introducing special treatments of the

boundaries, as detailed later.

ous field functions up to a certain order of accuracy.

They are related to the kernel moments and should

be satisfied both in the continuous and in the discrete

frameworks. More details can be found in l

*et alii*

*et alii*(2006).

motion are written in Lagrangian form, by using La-

grangian instead of Eulerian time derivatives.

proximation.

method.

**SPH NUMERICAL APPROACH IN MODELLING 2D MUDDY DEBRIS FLOW**

*Italian Journal of Engineering Geology and Environment - Book www.ijege.uniroma1.it © 2011 Casa Editrice Università La Sapienza*

speed is larger than the bulk velocity of flow is used,

then it is possible to keep density variations as low

as desired. The artificial equation of state used in this

paper is as follows:

*c*is the artificial speed of sound;

*ρ*

*0*

flow [see m

cantly affect the results, as long as the artificial speed of

sound is chosen to be large enough.

*VISCOSITY REGULARIZATION*

proaching to zero. It is therefore impossible to numer-

ically reproduce this behaviour in a straightforward

manner. In the calculations, the expression of eqn. (5)

for viscosity has been regularized and the infinite val-

ue that is obtained when the strain rate approaches to

zero has been replaced with a finite value. The expres-

sion used to regularize viscosity has been proposed by

P

*B*is related to the maximum vis-

strain rate is zero. In this case the maximum viscosity

value is given by:

*B*has been

cal behaviour of the fluid but small enough to avoid

prohibitively small time steps.

*BOUNDARY TREATMENT*

stress tensor divergence. Such an expression has been

proposed by C

ity too. It is possible to find a proof in Espanol et al.

(2003) that it accounts also for a bulk viscosity coeffi-

cient, which is equal to 5/3 times the dynamic shear vis-

cosity. Despite not having any bulk viscosity in the con-

stitutive law that is being used in the paper, we still think

it could be worth using the expression in eqn. (12) for

the deviatory tensor divergence for two main reasons:

if the velocity field is constant in space;

as the velocity field in the weakly compressible fluid

approximation is nearly divergence-free;

(10) and to use then the following SPH approximation

for its divergence:

actly reproduce the Herschel- Bulkley rheology within

the accuracy of the SPH approximation, without add-

ing any kind of bulk viscosity to the flow. Neverthe-

less, it lacks some good properties of the formulation

of eqn. (12) like the zero-th order consistency and the

angular momentum conservation

*ARTIFICIAL COMPRESSIBILITY*

which often leads to an increase of the computational

time. Therefore, it is more practical to approximate

the uncompressible medium with a weakly compress-

ible one, thus allowing the calculation of the pressure

from the density with a stiff equation of state which

introduces an artificial compressibility in the fluid.

*L. MINATTI & A. PASCULLI*

*5th International Conference on Debris-Flow Hazards Mitigation: Mechanics, Prediction and Assessment Padua, Italy - 14-17 June 2011*

errors introduced by neglecting the terms in the right

hand side of eqn. (9).

*et alii*(1996) used ghost

dition. Ghost particles have also been used in various

manners for particle approximations near boundaries

by t

*et alii*(1994), m

*et alii*(1997) and f

*et alii*(2009) by using point symmetry.

In this paper, boundaries have been treated by

is three times narrower than the fluid particles one.

Boundary particles prevent fluid particles from pass-

ing through the domain boundaries by exerting a nor-

mal force on them. They also interact with the fluid

particles via SPH summations through their viscosity,

thus enforcing the no-slip condition. The expression

used for the boundary forces has been proposed by

m

*et alii*(2009):

*k*

*b*

erted by the boundaries on the fluid;

Summation is extended to all the boundary particles lo-

cated within the support domain of fluid particle

*i*.

*et alii,*2009), that the

the direction parallel to the boundary, making boundary

forces being exerted only along the normal direction. It

is also a symmetrical formulation, thus conserving the

linear momentum of the particles system.

*TIME INTEGRATION*

*et alii*

serves the linear momentum of the particles system.

The time step

*Δt*is controlled by a C.F.L. condition

depending on the artificial speed of sound, the viscous

interactions between particles and on the interactions

with boundary particles, according to the following

equation:

The minimum time step value is sought over each

couple of interacting particle;

*h*

*ij*

*=(h*

*i*

*+h*

*j*

*)/2*is a symmetrised smoothing length be-

*c*

*ij*

*=(c*

*i*

*+c*

*j*

*)/2*is a symmetrised artificial speed of sound

*ρ*

*ij*

*(ρ*

*i*

*+ρ*

*j*

*)/2*is a symmetrised density between the pair

*dp*is the initial fluid particles spacing;

*β*is the ratio between the boundary and the fluid par-

ticles spacing;

The C.F.L. number, indicated as

*CFL*, has been set

equal to 0.5.

**MODEL TESTING**

experiments performed by l

*et alii*(1994) and

mudflow dam break problem in a laboratory flume, by

quickly opening a gate. The experimental setup they

used is briefly shown in the figure below:

sketched in the picture, recorded the mudflow wave

heights in time. The authors used waterclay mixtures

prepared in laboratory with different concentrations in

order to simulate mudflows. Herschel-Bulkley rheo-

logical parameters for the used mixtures were fitted

to measures, carried out with a rheometer by the au-

thors. They performed tests with four different kinds

of mixtures, named A, B, C and D. The rheological

*Fig. 1 - 3D sketch of the experimental setup in l*

*AiGle*

*et*

*alii (1994) and l*

*AiGle*

*(1997)*

**SPH NUMERICAL APPROACH IN MODELLING 2D MUDDY DEBRIS FLOW**

the number of particles required to fill the entire com-

putational domain gets higher. As the physics of the

flow is dominated by viscosity, the choice of the SPH

discretization for the momentum equation plays an im-

portant role. We have tried both eqn. (12) and (14) with

different particles spacing to reproduce the experiment

of Tab.1. Accuracy is improved at higher resolutions for

both equations but even though eqn. (12) seems to pro-

vide more accurate and less noisy results, it seems also

to be more dissipative than eqn. (14). The influence of

resolution and the issues regarding what would be the

most appropriate choice for the momentum equation

will be discussed further in a future paper

*dp*= 3 mm. The fol-

lowing figure shows the initial particle disposition,

right before the opening of the gate:

periments of l

*et alii.*(1994) and (1997) from

the thickness of the layer of particles passing below.

Finally, a comparison of our numerical wave heights

parameters are summarized in the following table:

were therefore completely determined by the height

of the material stored behind the gate, the flume slope

and the material characteristics. l

*et alii*(1994)

parameters controlling the resulting flow conditions.

The values of the scaling parameters in the tests they

performed can represent a wide range of real situa-

tions at smaller scales (C

of Tab. 1 represent a realistic natural material (with

*ρ*

*0*

*τ*

*c*

*k*= 290 Pa·s

*et alii,*1997). A field determina-

of a natural debris flow can be found in C

*et*

*alii*(1998), even though the values provided for the

parameters in the paper are in a slightly different range

than the ones indicated above. l

*et alii*(1994)

on the shallow water equations and solved it with a

Godunov conservative finite volumes scheme. We

tested our SPH model by simulating the experiment of

l

*et alii*(1997) using the rheological parameters

the experimental conditions indicated in Tab. 1.

placed to simulate the closed gate at the initial time

step. Particles were initially stored in a uniform lat-

tice according to their spacing and were given a

*ρ*

*0*

ter that, the gate opening was simulated by removing

the layer of boundary particles placed on it, thus re-

leasing the mudflow. The artificial Mach number was

set to 0.2, the boundary force constant of eqn. (18)

was set to

*k*

*b*

*= g*·

*H*

*g*= gravity acceleration,

*H*

*B*of eqn. (17) for the viscosity regularization

was set to 10 s.

*dp*, which plays the same role of the grid

spacing in finite differences schemes. A decrease in the

spacing improves the accuracy and reduces the numeri-

*Fig. 2 - Initial particles disposition (lengths are in meters*

*and particles are colour coded according to their*

*pressure in Pa)*

*Fig. 3 - Particles disposition after 0.15 s from gate open-*

*ing (lengths are in meters and particles are colour*

*coded according to the module of their velocity in*

*m/s)*

*L. MINATTI & A. PASCULLI*

*et alii*(1994) and

is capable of providing more information than the

one-dimensional model. Nevertheless, a preliminary

analysis shows that it is too diffusive, while still being

able to produce a good agreement with experimental

results. While the SPH model numerical diffusivity is

reduced by increasing the resolution, other kinds of

improvements haven’t been fully tested yet, but they

seem possible by using different SPH versions of the

momentum equation. In further papers, parametric

studies with comparison with actual debris flow phe-

nomena will be discussed.

*et alii*

numerical diffusion, as the calculated wave’s velocity

is lower than in the experimental data. This can be

improved either by further increasing the resolution or

by switching to a less diffusive SPH formulation for

the momentum equation.

**CONCLUDING REMARKS**

and has been tested to simulate some mudflow dam

*Fig. 4 - Particles disposition after 0.60 s from gate open-*

*ing (lengths are in meters and particles are colour*

*coded according to the module of their velocity in*

*m/s)*

*Fig. 5 - Time plot of the simulated wave heights at the lo-*

*cations of gage 1 (red line), gage 2 (blue line) and*

*gage 3 (black line). Experimental measurement s*

*of gages are represented by dots coloured as the*

*respective lines*

**REFERENCES**

*Continuum theories of mixtures: basic theory and historical development.*Q. J. Mech. Appl.

**29**(2)

**:**209-244.

*Computational Fluid Dynamics.*Cambridge University Press, New York.

*Modelling confined multi-material heat and mass flows using SPH*. Appl. Math. Modelling

**22:**981-993.

*Steady, laminar, flow of concentrated mud suspensions in open channel*. Journal of Hydraulic Research

**32**(4)

**:**535-559.

*Recognition, classification and mechanical description of debris flows*. Earth-Sci. Rev.

**40:**

*Direct determination of rheological characteristics of*

*debris flow*. Journal of Hydraulic Engineering 124(8): 865-868.

*Granular Flows and Numerical Mod-*

*eling of Landslides.*Debrisfall Assessment in Mountain Catchments for local end-users Contract N. EVG1-CT-1999-00007.

*Smoothed dissipative particle dynamics*. Physical Review

**67:**026705.

*A new 3D parallel scheme for free surface flows*. Computers & Fluids

**38:**1203-

*The physics of debris flows*. Reviews of Geophysics

**35**(3)

**:**245-296.

**227:**

**58:**

*The Particle Finite Element Method: a powerful tool to solve incompressible*

**SPH NUMERICAL APPROACH IN MODELLING 2D MUDDY DEBRIS FLOW**

*flows with freesurfaces and breaking waves*. Int. J. Numer. Meth. Engng.

**61:**964-989.

*Modélisation numérique des écoulements de laves torrentielles*. La Houille Blanche

**3:**50-56.

l

*SPH-based numerical investigation of mudflow and other complex fluid flow interac-*

*tions with structures*. Comput. Geosci. 11: 297-306.

*Smoothed Particle Hydrodynamics a meshfree particle method.*World Scientific Publishing.

*Restoring particle consistency in smoothed particle hydrodynamics*. Applied Numerical Math-

**56:**19-36.

*Simulating free surface flows with SPH.*Journal of Computational Physics

**110:**399-406.

*SPH and Riemann solvers.*Journal of Computational Physics

**136:**298-307.

*Smoothed particle hydrodynamics.*Reports on Progress in Physics

**68:**1703-1759.

*Modeling low Reynolds number incompressible flows using SPH.*Journal of Computa-

**136:**214-226.

*The Particle Finite Element Method*. An Overview. International

**1**(2)

**:**267-307.

*Flows of materials with yield.*Journal of Rheology

**31**(5)

**:**385.

*A rheologic classification of subaerial sediment-water flows*. Rev. Eng. Geol., Geol.Soc.

*Computational Granular Dynamics. Springer Berlin Heidelberg*, New York

*Smoothed particle hydrodynamics*: some recent improvements and applications. Com-

**138:**375-408.

*Debris Flow.*IAHR Monograph, published for IAHR by A.A.Balkema, Rotterdam

*Numerical simulation of viscous flow by smoothed particle hydrodynamics*.