# Polarization of contact forces in multi-contact systems

Serge Dumont, Jérôme Fortin, Youssef Ouafik

► To cite this version:

Serge Dumont, Jérôme Fortin, Youssef Ouafik. Polarization of contact forces in multi-contact systems. Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, 2010, 19 (1-3), pp.77-88. 10.3166/ejcm.19.77-88. hal-01113743

HAL Id: hal-01113743

https://hal.science/hal-01113743

Submitted on 12 Mar 2023

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Polarization of contact forces in
multi-contact systems

S. Dumont\* — J. Fortin\*\* — Y. Ouafik\*\*\*

\* LAMFA UPJV - CNRS UMR 6140

33 rue Saint-Leu, F-80000 Amiens, France
serge.dumont@u-picardie.fr

\*\* INSSET UPJV - LTI EA 3899

48 rue Raspail, F-02100 Saint-Quentin, France
jerome.fortin@u-picardie.fr

\*\*\* CEMEF - CNRS UMR 7635

BP 207, F-06904 Sofia Antipolis cedex, France
youssef.ouafik@mines-paristech.fr

RÉSUMÉ. Ce travail consiste à identifier des lois de comportement homogénéisées rendant compte du comportement global des systèmes multi-contacts. A l'heure actuelle, ces systèmes sont généralement analysés avec la Mécanique des Milieux Continus et la micro-mécanique ou approche multi-échelle. Cette étude a pour but de quantifier à l'aide de la méthode des Eléments Discrets (Jean, 1999; Fortin et al., 2005) les phénomènes de polarisation des efforts de contact.

ABSTRACT. The aim of this study is to identify the homogenized laws modeling the overall be-
havior of multi-contact systems. At the moment, these systems are generally analyzed either by
continuum mechanics or micro-mechanics and a multi-scale approach. These approaches dif-
fer from the phenomenological approach traditionally used for modeling the behavior of solid
materials which is based on mathematical formulations developed in the framework of thermo-
dynamics, whose constants are determined from results of laboratory tests. The lack of basic
physics in these formulations leads to mathematical models that are often complex and difficult
to identify. The multi-scale approach appears well suited to address these difficulties. This
study aims at quantification using the Discrete Element method (Jean, 1999; Fortin et al., 2005)
polarization phenomena of contact forces.

MOTS-CLÉS : Méthode des Eléments Discrets ; Méthode multi-échelle ; Frottement.

KEYWORDS: Discrete Element Method ; Multi-scale Method ; Friction.
# 1. Context

The aim of this study is to identify the homogenized behavior laws modeling the overall behavior of multi-contact systems. At the moment, these systems are generally analyzed either with continuum mechanics or micro-mechanics and multi-scale approach. These approaches differ from the phenomenological approach traditionally used for modeling the behavior of solid materials, which is based on mathematical formulations developed in the framework of thermodynamics, whose constants are determined from results of laboratory tests. The lack of basic physics in these formulations leads to mathematical models that are often complex and difficult to identify. The multi-scale approach appears well suited to address these difficulties. This study aims to quantify polarization phenomena of contact forces using the Discrete Element method (Jean, 1999; Fortin *et al.*, 2005). Polarization means that the main part of the contact forces in a relevant zone, for example on the walls of a silo, are on the boundary of the Coulomb cone.

# 2. Continuous approaches

Classically, to numerically study the behavior of a deformable body submitted to various solicitations, we use the Finite Element Method (FEM) which is an important tool for the analysis of structures. It is based on the Mechanics of Continuous Fields. Even though the bodies that are considered are not continuous, the hypothesis of continuity brings a simplification which makes possible the resolution of the problems of classic mechanics.

## 2.1. Stability of multi-contacts systems

Contrary to a liquid, a granular material can offer an oblique free surface. For a dry and non cohesive material, it is, however, impossible to tilt vertically : a pile of sand or earth cannot form a bank possessing a bigger slope than a certain critical angle according to the material and the geometry of the sand pile. Careful observation of the formation of a sand pile shows in fact that there are two critical angles : the angle of starting up $\theta_{start}$, and the angle of stopping $\theta_{stop}$. As soon as the angle of inclination exceeds $\theta_{start}$ an avalanche occurs : the material slides on the surface, which reduces the slope of the sand pile. This flow stops as soon as the angle of inclination becomes lower than $\theta_{stop}$.

The most direct interpretation would be to associate this coefficient of macroscopic friction with the microscopic friction between grains. Nevertheless, a pile of perfectly smooth grains, as soon as the first layer is maintained fixed, can also lead to an oblique free surface. The grain / grain friction is thus not necessary for the existence of a critical angle. Where then does this stability come from ? A first
element comes from the rigidity of the grains. To start moving, the pile has to dilate.

For dense flows on a free surface, the integrated approach of Saint-Venant is the most widely used. It consists in integrating the equations of preservation on the thickness of the layer into movement.

$$ \frac{\partial h}{\partial t}(x,t) = -\nabla \cdot q(x,t) + f(x,t) \quad [1] $$

This approach allows freedom the internal rheology of the layer in movement by considering only its averaged effects on the thickness of the layer in movement. The missing relations are generally chosen so as to report the studied phenomenon. We then obtain a complete set of equations which describes the dynamics of the granular material through the evolution of overall heights as the thickness of the layer in movement or the average flow. The flow of the granular material is connected in a thin boundary layer moving down slopes of the growing sand pile. There exists a function $m(x,t) \ge 0$ such that :

$$ q(x,t) = -m(x,t)\nabla h(x,t). \quad [2] $$

For numerical simulation, this equation has been discretized in time with an implicit Euler scheme, and a Finite Element Method is used for the space discretization of the domain Omega, Figure 1 (see (Dumont et al., 2009) for more details).

<figure><img src="image_2.png" /><figcaption>**Figure 1.** The sand source is localized upon the center of the domain</figcaption></figure>

However, this type of approach is purely descriptive and does not allow an understanding of the microscopic mechanisms responsible for continuous behaviors. On the scale of the grain, the two fundamental questions on the statics of granular assemblies concern the distribution of strengths when the assembly is at equilibrium and the stability conditions of this equilibrium.
## 2.2. Rheology of the ensiled materials

One of the important problems in the design of silos is the prediction of the vertical wall stress, induced by the stored granular material. This stress distribution depends of course on the granular material parameters and on the flow during the process of discharge. By modeling the ensiled granular medium as a continuous medium we can obtain relatively simple models for the calculation of the vertical wall stress. In 1885, Janssen (Janssen, 1895) offered a first model that provided a qualitative understanding of the saturation effect in granular silos. The idea of the Janssen model is that the vertical wall stress does not vary linearly with depth. The simplicity of the Janssen model (Janssen, 1895) comes because the granular media is considered as a continuous environment in a quasi-static state, that is that movements inside the silo are small enough for the global movement of the material to be considered as equal to zero. The main directions are in a perpendicular vertical plane on both walls of the silo (Figure (2)).

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Representation of a slice of ensiled matter</figcaption></figure>

The Janssen model also supposes the existence of a constant ratio $\lambda$ between the horizontal component of the stress tensor and its vertical component. This implies that $\lambda = \frac{\sigma_{xx}}{\sigma_{zz}}$ is constant in all the silo. On the other hand, the ensiled material is supposed to be in a state of break described by the Mohr-Coulomb criterion

$$ \tau_{xz} = \mu\sigma_{xx} \quad \text{with} \quad \mu = \tan \phi, \qquad [3] $$

where $\sigma_{xx}$, $\sigma_{zz}$ indicate the main stress and $\tau_{xz}$ the shear stress. Finally, the theory of Janssen is limited to the non cohesive, isotropic and homogeneous materials. Thus, Janssen was interested in the determination of the stress tensor assumed to depend only on $z$. The origin of the axis $z$ in the silo corresponds to the free surface of the material. On the other hand, Janssen considers that the upper surface of the ensiled material is without stress ($\sigma_{zz} = 0$ en $z = 0$). We finally obtain the stress expression
$$ 
\begin{aligned}
\sigma_{zz} &= \frac{\gamma R_h}{\mu \lambda} \left( 1 - \exp \left( -\frac{\mu \lambda}{R_h} z \right) \right) \\
\sigma_{xx} &= \lambda \sigma_{zz} \\
\tau_{xz} &= \mu \lambda \sigma_{zz}
\end{aligned}
\qquad [4]
$$

where $\gamma$ represents the volume weight, and $R_h$ indicates the hydraulic beam. However, the Janssen theory rests on certain contradictory hypotheses and contains some limitations; the Janssen theory assumes that $\sigma_{xx}$ and $\sigma_{zz}$ are the main stress connected by a constant coefficient in all the silo. These limitations of the Janssen theory explain the significant gaps existing between the stresses measured experimentally and those predicted by the Janssen theory. We notice that the hypothesis of continuity seems difficult to admit for systems which consist of several stiff or deformable parts, which are inter-connected. We then speak of multi-contact systems (Rahmoun, 2006).

### 3. Discrete approach

In fields where a collection of bodies between which the one-sided connections, usually affected by Coulomb friction, may become established or broken, these equations of constraints lead to a problem of non linear complementarity which cannot be solved by a linear programming method. The separation of surfaces in the case of sliding results from non compliance with the hypothesis of normality which implies a speed corresponding to the dilation of the interface.

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Multi-contact systems</figcaption></figure>

To each couple of particles $\Omega_i$ and $\Omega_j$ which may enter in contact, is associated a local reference whose axes are oriented according to the two unit vectors $\mathbf{n}$ and $\mathbf{t}$, respectively, normal and tangential vectors in the contact plane. The normal $\mathbf{n}$ is
directed from $\Omega_j$ to $\Omega_i$. The dual variables are $\dot{\mathbf{u}}$, the local relative velocity between $\Omega_i$ and $\Omega_j$, and the contact reaction $\mathbf{r}$ of $\Omega_j$ on $\Omega_i$. In the local base, they are written by

$$ \dot{\mathbf{u}} = \dot{\mathbf{u}}_t + \dot{u}_n \cdot \mathbf{n}, \quad \mathbf{r} = \mathbf{r}_t + r_n \cdot \mathbf{n}, \qquad [5] $$

where $\dot{u}_n$ is the normal velocity, $\dot{\mathbf{u}}_t$ the sliding velocity, $r_n$ the contact pressure and $\mathbf{r}_t$ the friction force. The introduction of Coulomb friction $K_\mu$ defined by

$$ K_\mu = \{(r_n, \mathbf{r}_t) \text{ such that } f(r_n, \mathbf{r}_t) = \|\mathbf{r}_t\| - \mu r_n \le 0\} \qquad [6] $$

leads to a non-linear problem which cannot be solved by a linear programming method. It is possible to prove that the Coulomb unilateral contact law with dry friction does not satisfy the cyclic monotonic condition of Rockaffelar. Thus there does not exist an associated formulation in terms of sub-differential using a pseudo-potential. However one can associate a bi-potential as follows :

$$ b_c(-\dot{\mathbf{u}}, \mathbf{r}) = \Psi_{\mathbb{R}^{-}}(-\dot{u}_n) + \Psi_{K_\mu}(\mathbf{r}) + \mu r_n \| - \dot{\mathbf{u}}_t \|. \qquad [7] $$

The condition of non inter-penetrability $\dot{u}_n \ge 0$ is represented by the indicatory function of $\mathbb{R}^{-}$, noted $\Psi_{\mathbb{R}^{-}}(-\dot{u}_n)$, which is equal to zero when $-\dot{u}_n \le 0$ and to $+\infty$ otherwise. The contact bi-potential also takes infinite values if the condition $\mathbf{r} \in K_\mu$ is not satisfied. This bi-potential of contact is bi-convex (convex with respect to each of the variables) and satisfies :

$$ \forall -\dot{\mathbf{u}}, \mathbf{r} \in \mathbb{R}^3, \quad b_c(-\dot{\mathbf{u}}, \mathbf{r}) \ge -\dot{\mathbf{u}} \cdot \mathbf{r}. \qquad [8] $$

Moreover, the couples for whom equality is reached in the previous relation, are called extremal couples :

$$ b_c(-\dot{\mathbf{u}}, \mathbf{r}) = -\dot{\mathbf{u}} \cdot \mathbf{r} \Leftrightarrow \mu r_n \| -\dot{\mathbf{u}}_t \| = -(\dot{\mathbf{u}}_t \cdot \mathbf{r}_t + \dot{u}_n r_n). \qquad [9] $$

These couples verify the Coulomb unilateral contact law with dry friction and the inverse law, which can be written implicitly :

$$ -\dot{\mathbf{u}} \in \partial_r b_c(-\dot{\mathbf{u}}, \mathbf{r}), \quad \mathbf{r} \in \partial_{-\dot{\mathbf{u}}} b_c(-\dot{\mathbf{u}}, \mathbf{r}) \qquad [10] $$

where $\partial_x b_c$ denotes the sub-differential of $b_c$ with respect to the variable $x$. Classically the resolution of the Coulomb unilateral contact law with dry friction needs two
principles of minimization : the first one for unilateral contact, and the second one for friction. The use of the bi-potential only needs one variational principle where the contact and the friction are coupled. For the resolution of the unilateral contact law with bi-potential formalism, we use the augmented Lagrangian method. First, let us write a relation as follows :

$$ \forall \mathbf{r}' \in K_{\mu}, \quad b_c(-\dot{\mathbf{u}}, \mathbf{r}') - b_c(-\dot{\mathbf{u}}, \mathbf{r}) \geq -\dot{\mathbf{u}}(\mathbf{r}' - \mathbf{r}). \quad [11] $$

Now let us choose a positive arbitrary coefficient $\rho$, whose value will be fixed later to ensure the numerical convergence of the algorithm. Then inequality (11) can be written :

$$ \forall \mathbf{r}' \in K_{\mu}, \quad \rho b_c(-\dot{\mathbf{u}}, \mathbf{r}') - \rho b_c(-\dot{\mathbf{u}}, \mathbf{r}) + [\mathbf{r} - (\mathbf{r} + \rho(-\dot{\mathbf{u}}))].(\mathbf{r}' - \mathbf{r}) \geq 0. \quad [12] $$

Using now the definition (7) of the contact bi-potential, relation (12) becomes with $\dot{u}_n \geq 0$ and $\mathbf{r} \in K_\mu$ :

$$ \forall \mathbf{r}' \in K_{\mu}, \quad (\mathbf{r} - \tau).(\mathbf{r}' - \mathbf{r}) \geq 0 \quad [13] $$

where

$$ \tau = \mathbf{r} - \rho[\dot{\mathbf{u}}_t + (\dot{u}_n + \mu|| - \dot{u}_t||).\mathbf{n}] $$

denotes the augmented reaction. Relation (13) implies that $\mathbf{r}$ is the projection of $\tau$ onto Coulomb cone :

$$ \mathbf{r} = \mathrm{proj}(\tau, K_{\mu}) \quad [14] $$

and can be solved with a Usawa-like algorithm; let $(-\dot{\mathbf{u}}^i, \mathbf{r}^i)$ be an approximation of $(-\dot{\mathbf{u}}, \mathbf{r})$ at the iteration i. Then the calculus of $\mathbf{r}^{i+1}$ is split into two steps :

$$ \text{prediction :} \quad \tau^{i+1} = \mathbf{r}^i - \rho[\dot{\mathbf{u}}_t^i + (\dot{\mathbf{u}}_n^i + \mu|| - \dot{\mathbf{u}}_t^i||).\mathbf{n}], \quad [15] $$

$$ \text{correction :} \quad \mathbf{r}^{i+1} = \mathrm{proj}(\tau^{i+1}, K_{\mu}). $$

where the projection of Coulomb's cone leads, according to the value of $\tau$, to one of these states : no contact, contact with adherence or sliding contact, (Figure 4). Conventionally, at each time step, the contact forces in the system are determined repeatedly by the method of successive balances based on a Gauss-Seidel algorithm for the 2D version. Each contact force is calculated by adopting temporary values over the other contacts. Convergence is obtained when the force confirms the unilateral contact law with dry friction.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** *Projection of the augmented reaction onto the friction cone*</figcaption></figure>

**4. Homogenized laws of behavior**

While the stress tensor has for a long time been defined for a continuous environ-
ment, the definition of an average stress tensor from the strengths being applied on the
scale of the grain itself is much less clear for granular media. On one hand, there are
several different definitions which do not generally lead to a symmetric average stress
tensor. We can find a comparison between the various expressions obtained as well as
an analysis of an order of height of the antisymmetric part of the average stress tensors
defined. Eventually, it is important to note that all these works are based on a static
analysis of the problem, and are valid only for a quasi-static arrangement of granular
media.

<figure><img src="image_6.png" /><figcaption>**Figure 5.** *Homogenization procedures : to define a Representation Elementary Vo-lume*</figcaption></figure>
In the dynamic state, only J.J. Moreau proposes a generalization of the previous definitions based on the notion of internal moment. Hence, two problems still remain open ; the first concerns the choice of a substantial definition of the average stress tensor in a granular medium in quasi-static conditions. The second concerns an extension of this definition in dynamics. To try to clarify the problem, we have proposed a general definition of the average stress tensor for a granular medium, valid in statics as well as in dynamics. This definition, which follows (de Saxcé et al., 2004) takes into account not only the strengths of contact, but also the strengths of volume (gravity and slowness) being applied on the scale of grains. The average stress tensor so defined possesses the properties of a Cauchy stress tensor : it is symmetric and invariant by translation.

$$ \Sigma_p = \frac{1}{V_b} \left( \int_{V_b} \vec{x} \otimes \left[ \rho(\vec{g} - \ddot{\vec{x}}) \right] dV + \vec{x}_c \otimes \vec{r} \right) = \Sigma_r + \Sigma_g + \Sigma_\gamma \quad [16] $$

where $\Sigma_r$, $\Sigma_g$ and $\Sigma_\gamma$ denote the respective contributions of the reaction force, of the gravity force and of the acceleration forces. We obtain, after computations, (de Saxcé et al., 2004).

$$ \Sigma_p = \rho g a \left\{ \left( \begin{array}{cc} 0 & \frac{\sin\varphi}{6} \\ \frac{\sin\varphi}{6} & -\cos\varphi \end{array} \right) + \frac{\mathcal{L}}{a} \left( \begin{array}{cc} \frac{\sin\varphi}{3} & 0 \\ 0 & \frac{\sin\varphi}{3} \end{array} \right) \right\} \quad [17] $$

where $\mathcal{L}$ can be interpreted as the mean free path of the bead between two collisions.

# 5. Application

## 5.1. Test case

Quasi-static examples, presented in Fig 6, are considered. A regular array of equal-size disks confined in a box is considered. The box is composed of two vertical and two horizontal flat walls. The particles as well as the walls are assumed to be perfectly rigid. Each particle is subjected to the gravitation force $g$ and to contact force resulting from neighborhood particles and walls. For three different geometries, the normal contact forces for each particle are computed. The assumptions and results are displayed in Figure 6.

The widths of inter-center segments are proportional to the corresponding normal contact force intensity. For the rectangular arrangement, for obvious physical reasons, the normal force increases when the height of the disk center with respect to the base of
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Normal contact forces(a) rectangular - (b) triangular - (c) cannon ball</figcaption></figure>

the box decreases. For the triangular arrangement, we observe that contacts between the particles and the vertical wall exist. This situation is typical regarding the complex behavior of granular materials. By comparison with the analytical computation, we obtain the same values concerning the normal contact forces. For the cannon-ball arrangement, friction, which is not displayed here, is essential to the stability of the system.

## 5.2. Polarization of contact forces

First, we consider a silo where the properties of the contact stresses on the lateral walls allow an explanation of the macroscopic properties, such as the stopping of the flow at the bottom of the silo, or the modeling of the stresses in the silo with the Janssen model. In the experiment presented below, we have tried to determine if the contact stresses on the vertical walls are in the so-called Coulomb conditions, i.e. the normal and the tangential reaction are on the slip surface of the Coulomb cone. We have considered here a sample of 350 particles, with a radius equal to 1.5mm (±10%) in a silo. The static friction coefficient is equal to 0.25. The first part of the experiment was to settle the sample under the influence of gravity, and then to relax the system by lowering the bottom wall in order to put the system in the Coulomb condition (Figure 7a). At the end of the relaxation, we compute the average constraint on the lateral walls, and compare the results with those given by the Janssen approach (Figure 7b).

In a second step, we have considered a shear flow, Figure 8, where the quantification of the polarization will permit us to propose a model of boundary layer for this type of flow, like what has been achieved in the case of masonry (Lebon *et al.*, 2008). This type of study has already been done with the DEM in the case of non-sheared media, as for example in a rotating drum (Renouf *et al.*, 2005).

In the first term of (17), the $\sigma_{yy}$ component is negative, which corresponds, with our sign convention, to a compression stress. It can be interpreted as due to the contact reaction onto the bead (induced by the action of gravity) which prevents it from penetrating through the plane. The $\sigma_{xy}$ component of the first term represents the friction effect. In the second term of (17), the resulting hydrostatic stress is positive, which
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Polarization of contact forces - (b) Comparison of Janssen - DEM</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Shear flow</figcaption></figure>

corresponds to a traction state due to the centrifugal effects. We call it the centrifugal stress tensor. However, in a dense granular medium, because of frequent collisions, the mean free path of particles is small with respect to their size. Hence, the ratio $\underline{L}/a$ remains small with respect to 1 and the second term in (17) can be neglected. It is worthwhile noting that both contact reaction and gravity produce non-symmetrical tensors $\Sigma_r$, $\Sigma_g$. This shows that it is important not to neglect the inertia tensor $\Sigma_I$ in the calculation of the mean stress tensor. As the inertia forces have to balance the other forces, contact reactions and gravity, they have the same order of magnitude.
## 6. Conclusion

In this study we have shown both for a quasi-static example as in the silo and for a dynamic example as in a shear flow over a sandpile that the position of the contact forces on the Coulomb cone enables to explain the overall behavior of the granular material.

Aknowledgement : this program is supported by the French program ANR JCJC 41831.Thanks to Claude James for reviewing english of the paper.

## 7. Bibliographie

Chetouane B., Approche combinée d'éléments finis/éléments discrets pour la modélisation des structures maçonnées, Thèse de doctorat, Université de Montpellier 2, 2004.

de Saxcé G., Fortin J., Millet O., « About the numerical simulation of the dynamics of granular media and the definition of the mean stress tensor », *Mechanics of Materials*, vol. 36, p. 1175-1184, 2004.

de Saxcé G., Feng Z. Q., « New inequality and functional for contact with friction : The implicit standard material approach », *Mechanics of Structures and Machines*, vol. 19, no 3, p. 301-325, 1991.

Dumont S., Igbida N., « On a Dual Formulation for the Growing Sandpile Problem », *European Journal of Applied Mathematics*, vol. 20, p. 169-185, 2009.

Fortin J., Simulation numérique de la dynamique des systèmes multicorps appliquée aux milieux granulaires, Thèse de doctorat, Université de Lille 1, 2000.

Fortin J., Millet O., de Saxcé G., « Numerical simulation of granular materials by an improved discrete element method », *Int. J. Numer. Meth. Engng*, vol. 62, p. 639-663, 2005.

Janssen H. A., « Versuche übergetreichedruck in silozellen », *Zeitschriftverein deutscher ingenieur*, vol. 39, p. 1045-1049, 1895.

Jean M., « The Non Smooth Contact Dynamics Method », *Compt. M. Appl. Math. Engrg*, vol. 177, p. 235-257, 1999.

Lebon F., Rizzoni R., « Asymptotic study of soft thin layer : the non convex case », *Mechanics of Advanced Materials and Structures*, vol. 15, p. 12-20, 2008.

Moreau J. J., « Sur les lois de frottement, de plasticité et de viscosité », *C. R. Acad. Sci. Série 1*, vol. 271, p. 608-611, 1970.

Rahmoun J., Etude théorique et numérique du comportement des milieux granulaires, Thèse de doctorat, Université de Lille1, 2006.

Renouf M., Bonamy D., Dubois F., Alart P., « Numerical simulation of two-dimensional steady granular flows in rotating drums : on surface flow rheology », *Physics of fluids*, vol. 17, p. 13303-13315, 2005.