# Modeling the impact of mudflows on protection structures

Matthias Preisig, Thomas Zimmermann

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

To cite this version:

Matthias Preisig, Thomas Zimmermann. Modeling the impact of mudflows on protection structures. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01422244

HAL Id: hal-01422244

https://hal.science/hal-01422244

Submitted on 24 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Modeling the impact of mudflows on protection structures

M. Preisig¹, Th. Zimmermann¹

¹ Laboratoire de simulation en mécanique des solides, EPFL
Station 18, CH-1015 Lausanne
matthias.preisig@epfl.ch, zimmermann@zace.com

**Résumé** — Un modèle biphasé de l’écoulement d’une mixture dans un repère mobile est proposé. Après chaque mise à jour des nœuds des deux phases un nouveau maillage est généré. Le calcul des fractions volumiques des deux phases est découplé du calcul des vitesses et des pressions. La description de l’écoulement comme un continu permet d’effectuer l’analyse d’intéraction fluide-structure. Sur le problème de l’impacte sur un obstacle d’une coulée de boue nous mettons en évidence la performance de la méthode pour le calcul de forces agissantes sur des structures.

**Mots clés** — Interaction fluide-structure, dangers naturels, écoulement biphasé.

## 1 Introduction

Every year throughout the world, debris flows, mudflows and landslides cause an immense amount of damage to property and people. Mitigation of these risks is two-fold : One, mitigation requires better zoning and urban planning measures to single out danger zones. Two, it requires a better understanding of these events in order to protect the existing infrastructure from damage. In both cases, numerical modeling can play an important role in risk mitigation by one, helping to establish maps that outline areas which are in the flow path, and two, computing characteristics of the event such as depth of the flow or slide and forces, both necessary for designing effective barriers, walls, nets or other protection devices.

Mudflows are flows of a mixture of water with soil, rocks and other debris. From a modeling point of view, they are characterized as gravity-driven flows of a multiphase material. During the event, the material undergoes large motion, which causes the boundary to change continuously. The boundary, and in particular the free surface, is not known a priori, but part of the solution.

Current state-of-the-art numerical models of mud- and debris flow use the thin-layer assumption, stating that the extension of the flow in and perpendicular to the direction of flow is much larger than the depth of the flow. This justifies the application of depth-averaging, where the flow is projected onto a line for simple channel models or onto a two-dimensional mesh (a map). The depth of the flow and the velocity are part of the solution in each mesh node. A hypothesis on the vertical distribution of stresses is required to integrate the governing equations over depth.

Depth-averaged models have proven to be very valuable for evaluating run-out distances and flow paths in order to edit hazard maps. Due to the thin-layer assumption, however, stresses cannot be predicted with very high accuracy and the application of these models for obtaining loads acting on protection devices is very limited. In order to obtain results that have the required degree of reliability, the next generation numerical model has to be full scale three dimensional. The present work is an important step in this direction.

The goal of this work is to propose a model for mud- or debris flows that allow the simulation of its impact on an obstacle, and to obtain detailed time histories of forces acting on the struc-
ture. The model has to be able to take into account the complex material behavior of a soil-water mixture. The algorithmic framework for describing the motion of the mixture has to allow the inclusion of structures in order to, in a next step, perform fully coupled fluid-structure interaction analysis.

Full scale modeling of a fluid undergoing large motion and whose boundary changes continuously poses special requirements to the spatial description. We chose a formulation where particles move according to their mass and acceleration while being subject to forces from the surrounding fluid. The incompressible Navier-Stokes equations are solved in a Lagrangian reference frame at each time step. The nodes are updated (moved) to the current spatial position after each step.

It is well known that Lagrangian methods suffer from problems related to mesh distortion. This can be dealt with by re-zoning the nodes and reconnecting the mesh. However, in the framework of simulating free-surface flows it is our opinion that the advantages of Lagrangian meshes outweigh their disadvantages.

The main difficulty in developing an updated Lagrangian method for simulating a two-phase mixture consists in finding an appropriate way to update the spatial coordinates of the nodes. Since in each node the velocities of two phases are computed, the update produces two new nodes, one for each phase. Multiplying the number of nodes by two at each Lagrangian update is not a solution for obvious reasons. We deal with this fundamental problem by re-creating a new set of nodes, on which the solution transported by the nodes at their updated spatial coordinates is interpolated. This approach solves at the same time the problem of mesh distortion.

## 2 Two-phase mixture

The material behavior of a mudflow is governed by the properties of its constituents, water and granular material. The grain sizes of the granular material range from very fine particles, as in clay, to gravel and even boulders. We restrict ourselves in this work to material with an important fine content, a characteristic property of sediments deposited by the melting of a glacier. The high fine content gives the mixture a muddy, viscous behavior. Keeping in mind that the algorithmic framework is not specific to any particular kind of material, the scope of this work is limited to a two-phase material where both the fluid and the solid phase are modeled as viscous fluids.

Both phases are smeared over the fluid domain using mixture theory. This means that no phase interfaces are considered. The presence of either phase is specified by its volume fraction. The volume fraction of a phase $p$ is given by

$$C_p = \frac{V_p}{V} \qquad (1)$$

where $V$ is a control volume and $V_p$ is the volume within $V$ that is occupied by phase $p$. $p$ is a placeholder which can stand for the solid phase ($p = s$) or the fluid phase ($p = f$). We assume that the control volume is entirely occupied by the two phases, therefore $C_s + C_f = 1$ holds. Considering a control volume that is fixed in space another important observation, which we will use later, can be made :

$$\frac{\partial C_s}{\partial t} + \frac{\partial C_f}{\partial t} = 0 \qquad (2)$$

## 3 Governing equations of two-phase flow

The governing equations of two-phase flow are one equation of total mass conservation and two equations of momentum conservation, one for each phase. The pressure $p$ is the same in both
phases, which is a reasonable assumption for mixtures with uniform grain size distributions. The stress in each phase is given by a constitutive relation and interaction between the two phases is governed by a momentum exchange relation. The model of the two phases smeared over the entire volume occupied by the mixture requires volume averaging of the governing equations of each phase. We follow the procedure of volume averaging as it is laid out in the book by Soo [1]. The detailed development can be found in Preisig [2]. The equation of conservation of momentum of a phase p is

$$C_p \rho_p \frac{D\mathbf{v}_p}{Dt} = \nabla \cdot (\sigma_p) + C_p \rho_p \mathbf{b} + K'_{\text{drag}} (\mathbf{v}_p - \mathbf{v}_q) \quad (3)$$

$\rho_p$ is the (intrinsic) density of phase p and $\mathbf{v}_p$ its (material) velocity. $\mathbf{v}_q$ represents the velocity of the other phase. $\langle \sigma_p \rangle$ is the stress in phase p, averaged over a control volume. The last term represents the momentum exchange between phases p and q. The drag force coefficient $K'_{\text{drag}}$ depends on the solid volume fraction $C_s$ and the viscosity of the fluid phase $\mu_f$: $K'_{\text{drag}} = K_{\text{drag}} \frac{C_s}{(1-C_s)^2} \mu_f$. $K_{\text{drag}}$ is a material parameter depending on the fluid and the solid phase.

The constitutive relations of the two-phase fluid are derived from the constitutive relation of a single-phase viscous fluid. The averaged total stress in the mixture is the sum of the averaged stresses in the individual phases:

$$\langle \sigma \rangle = \langle \sigma_s \rangle + \langle \sigma_f \rangle \quad (4)$$

This relation is satisfied for the following constitutive equations:

$$\langle \sigma_s \rangle = C_s (\tau(\mathbf{v}_s) + p\mathbf{I}) \quad \text{and} \quad \langle \sigma_f \rangle = C_f (\tau(\mathbf{v}_f) + p\mathbf{I}) \quad (5)$$

The deviatoric stress tensors are defined as

$$\tau(\mathbf{v}_s) = 2\mu_s \left( \dot{\mathbf{\epsilon}}(\mathbf{v}_s) - \frac{1}{3}(\nabla \cdot \mathbf{v}_s)\mathbf{I} \right) \quad (6)$$

$$\tau(\mathbf{v}_f) = 2\mu_f \left( \dot{\mathbf{\epsilon}}(\mathbf{v}_f) - \frac{1}{3}(\nabla \cdot \mathbf{v}_f)\mathbf{I} \right) \quad (7)$$

where $\mu_s$ and $\mu_f$ are the dynamic viscosities of the solid and the fluid phase. The rates of deformation are given by $\dot{\mathbf{\epsilon}}_s = \frac{1}{2}(\nabla \mathbf{v}_s + (\nabla \mathbf{v}_s)^T)$ and $\dot{\mathbf{\epsilon}}_f = \frac{1}{2}(\nabla \mathbf{v}_f + (\nabla \mathbf{v}_f)^T)$.

The equations of conservation of mass, written in an Eulerian reference frame and for incompressible behavior of both phases, are

$$\frac{\partial C_s}{\partial t} + \nabla \cdot (C_s \mathbf{v}_s) = 0 \quad \text{and} \quad \frac{\partial C_f}{\partial t} + \nabla \cdot (C_f \mathbf{v}_f) = 0 \quad (8)$$

The equation of conservation of mass of the mixture is obtained by combining the equations of conservation of mass of each phase and making use of Equation 2:

$$\nabla \cdot (C_s \mathbf{v}_s) + \nabla \cdot (C_f \mathbf{v}_f) = 0 \quad (9)$$

The boundary value problem consists in solving the following equations for velocities and
pressure, given the boundary conditions $\mathbf{g}_s$, $\mathbf{g}_f$, $\mathbf{h}_s$, $\mathbf{h}_f$ and the initial conditions $\mathbf{v}_{s,0}$, $\mathbf{v}_{f,0}$ and $p_0$:

$$C_s \rho_s \frac{D\mathbf{v}_s}{Dt} = \nabla \cdot [C_s(\tau(\mathbf{v}_s) + p\mathbf{I})] + C_s \rho_s \mathbf{b} + \mathbf{m}_{sf} \quad \text{on } \Omega \times ]0, T[ \quad (10)$$

$$C_f \rho_f \frac{D\mathbf{v}_f}{Dt} = \nabla \cdot [C_f(\tau(\mathbf{v}_f) + p\mathbf{I})] + C_f \rho_f \mathbf{b} + \mathbf{m}_{fs} \quad \text{on } \Omega \times ]0, T[ \quad (11)$$

$$0 = \nabla \cdot (C_s \mathbf{v}_s) + \nabla \cdot (C_f \mathbf{v}_f) \quad \text{on } \Omega \times ]0, T[ \quad (12)$$

$$\mathbf{v}_s = \mathbf{g}_s \quad \text{on } \partial\Omega_{g_s} \times ]0, T[ \quad (13)$$

$$\mathbf{v}_f = \mathbf{g}_f \quad \text{on } \partial\Omega_{g_f} \times ]0, T[ \quad (14)$$

$$\boldsymbol{\sigma}_s \cdot \mathbf{n} = \mathbf{h}_s \quad \text{on } \partial\Omega_{h_s} \times ]0, T[ \quad (15)$$

$$\boldsymbol{\sigma}_f \cdot \mathbf{n} = \mathbf{h}_f \quad \text{on } \partial\Omega_{h_f} \times ]0, T[ \quad (16)$$

$$\mathbf{v}_s(t=0) = \mathbf{v}_{s,0} \quad \text{on } \Omega \quad (17)$$

$$\mathbf{v}_f(t=0) = \mathbf{v}_{f,0} \quad \text{on } \Omega \quad (18)$$

$\partial\Omega_{g_p}$ denotes the part of the boundary on which we impose the displacement $\mathbf{g}_p$, while $\partial\Omega_{h_p}$ denotes the Neumann part, where we impose surface tractions $\mathbf{h}_p$.

# 4 Spatial discretization

The weak form of the above boundary value problem is discretized in space using linear triangular finite elements for both the velocity and the pressure. In order to avoid pressure oscillations caused by non-satisfaction of the LBB condition the weak form has to be stabilized by adding stabilization terms. Details can be found in [2].

# 5 Time integration scheme

Time stepping is performed using a generalized trapezoidal algorithm. Discretization in time of the semidiscrete matrix form at $t = t_{n+1}$ leads to:

$$\mathbf{M}(\mathbf{x}_{n+1}, C_s^{n+1}, C_f^{n+1})\mathbf{a}_{n+1} + \mathbf{K}(\mathbf{x}_{n+1}, C_s^{n+1}, C_f^{n+1})\mathbf{v}_{n+1} = \mathbf{F}_{n+1}^{\text{ext}} \quad (19)$$

$\mathbf{M}$ denotes the mass matrix, $\mathbf{K}$ the stiffness matrix, $\mathbf{a}$ the (material) acceleration and $\mathbf{F}^{\text{ext}}$ the external force vector. The algorithm consists of a predictor step, a step where the above linear system is solved for velocities and pressure, and a corrector step. At the end of both the predictor and the corrector step the nodes are updated to their new positions, the mesh is updated, the variables are mapped from the previous to the new mesh and the volume fractions $C_s$ and $C_f$ are recomputed.

The Lagrangian update of the nodes requires special attention due to the presence of two phase velocities in each node. The update of one node yields two new nodes, one for each phase. In order to avoid an exponential increase of the number of nodes a regularly distributed set of nodes is created, new elements are created using a Delaunay triangulation, and the variables are mapped onto this new mesh by linear interpolation. The nodes on the boundary are updated by a displacement increment that is computed from the velocity of the mixture. This assures satisfaction of Equation 9.

# 6 Computation of volume fractions

New volume fractions $C_s$ and $C_f$ are evaluated after each mesh update. While global volume conservation is assured by Equation 9 the conservation of the mass of each individual phase has to be taken care of when the new volume fractions are computed.
In a few words, the method can be summarized as follows : the Lagrangian update of the nodes transports the masses of the solid and the fluid phase to their new position. At the nodes of the new (regenerated) mesh the volume fractions can be computed by evaluating the solid and the fluid masses at the updated nodal positions of the previous mesh. This is done by an approximation of the local density of solid and fluid mass. We approximate this density by comparing the relative movement of surrounding nodes of the previous mesh with respect to a node of the new mesh.

# 7 Numerical results

In the following we show results of a simulation of the downhill propagation of a two-phase mixture and its impact on an obstacle. The obstacle, representing a protection dam, is modeled as a solid block which is placed at the bottom of a slope. The flow is initiated by the sudden release of a homogeneous two-phase mixture on top of a slope. Material parameters together with details of the discretization are given in Table 1. For comparison the same problem is simulated with a single-phase fluid, using the average material properties of the two-phase mixture.

<table><caption>Table 1 – Parameters of the simulation of a mudflow impacting on an obstacle.</caption><thead><tr><th></th><th>ρ<sub>s</sub></th><th>ρ<sub>f</sub></th><th>μ<sub>s</sub></th><th>μ<sub>f</sub></th><th>K<sub>drag</sub></th><th>C<sup>init</sup><sub>s</sub></th><th>Δt</th><th>No. of nodes</th></tr></thead><tbody><tr><td>Two-phase model</td><td>1000</td><td>500</td><td>100</td><td>2</td><td>10'000</td><td>0.5</td><td>0.002</td><td>≈ 1550</td></tr><tr><td>Single-phase model</td><td>750</td><td></td><td>51</td><td></td><td>-</td><td>-</td><td>0.002</td><td>≈ 1550</td></tr></tbody></table>

The outline of the two-phase fluid during the event is shown in Figure 1. Colors indicate solid volume fractions. The solid phase initially accumulates at the base, while the fluid phase stays

<figure><img src="image_3.png" /><figcaption>Figure 1 – Solid volume fractions on a mudflow impacting on an obstacle.</figcaption></figure>
on the surface and accumulates at the front of the flow. After the flow tip reaches the obstacle the solid phase quickly catches up with the faster flowing fluid phase, filling the space behind the barrier. The tip of the mixture shooting over the barrier is essentially fluid, due to lower viscosity and density. At the end of the simulation, that is after 1000 time steps, the total volume of the mixture has increased by 2.3%. We consider this error very small, given the relatively simple contact algorithm used in the model.

From the simulation the resultant force acting on an obstacle that obstructs the flow path is extracted. In Figure 2 the resultant force of the single-phase and the two-phase model are compared. The force is computed by integrating the pressure along the front side of the obstacle. Right after

<figure><img src="image_4.png" /><figcaption>Figure 2 – Resultant force acting on obstacle.</figcaption></figure>

the impact the force attains its peak, before it slowly decays to the hydrostatic level. The peak force right after the impact is higher in the case of a single-phase material. This can be explained with the lower density of the fluid phase, which reaches the obstacle first in the two-phase simulation. After a while however the difference vanishes. The fluid phase acts as a buffer, attenuating the impact of the denser solid phase.

## 8 Conclusions

A new and innovative numerical method for simulating two-phase free-surface flows has been developed in this work. The method is capable of simulating a wide range of problem types from sedimentation of solid material in a fluid (see [2]) to gravity-driven free-surface flows. The key feature of the method is an algorithm that allows the motion of two different constituents of a mixture to be followed in a Lagrangian reference frame.

In contrast to existing debris- or mudflow models the method implements a continuum approach, which allows to obtain detailed time histories and profiles of stresses, velocities or volume fractions. The algorithmic framework is kept general in order to allow any kind of constitutive model to be included. We expect the method to find a wide range of applications not only in the field of geophysical flows, but in any kind of problem involving the motion of two phases where interaction between the phases cannot be neglected.

## References

[1] S. L. Soo. *Multiphase Fluid Dynamics*. Science Press, Beijing, China, 1990.

[2] M. Preisig. *Modeling two-phase flows on moving domains*. PhD thesis, Ecole Polytechnique Fédérale de Lausanne, Lausanne, Switzerland, 2008. Thesis no. 4250.