Miguel A Fernández
email: miguel.fernandez@inria.fr

Mikel Landajuela
email: mikel.landajuelalarma@inria.fr

Numerical Analysis/Mathematical Problems in Mechanics Splitting schemes for incompressible fluid/thin-walled structure interaction with unfitted meshes 1

Splitting schemes for incompressible fluid/thin-walled structure interaction with unfitted meshes

Miguel Angel Fernández, Mikel Landajuela

ont été réalisés dans le développement et l'analyse de schémas de couplage explicites (à savoir, qui ne résolvent le fluide et le solide qu'une seule fois par pas de temps) pour l'approximation de (1)- [START_REF] Bukač | Fluid-structure interaction in blood flow capturing non-zero longitudinal structure displacement[END_REF]. Tous ces travaux (see, e.g., [START_REF] Burman | Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility[END_REF][START_REF] Guidoboni | Stable loosely-coupled-type algorithm for fluid-structure interaction in blood flow[END_REF][START_REF] Fernández | Incremental displacement-correction schemes for incompressible fluid-structure interaction: stability and convergence analysis[END_REF][START_REF] Fernández | Explicit Robin-Neumann schemes for the coupling of incompressible fluids with thin-walled structures[END_REF][START_REF] Burman | Explicit strategies for incompressible fluid-structure interaction problems: Nitsche type mortaring versus Robin-Robin coupling[END_REF][START_REF] Bukač | Fluid-structure interaction in blood flow capturing non-zero longitudinal structure displacement[END_REF][START_REF] Lukáčová | Kinematic splitting algorithm for fluid-structure interaction in hemodynamics[END_REF]) considèrent des maillages compatibles (Figure 1, gauche). Or, il est bien connu que pour beaucoup d'applications cette hypothèse devient critique très rapidement (grandes déviations de l'interface, structures en contact, etc). Dans le cadre de maillages non compatibles (Figure 1, droite) des schémas de couplage explicite ont été proposés et analysés dans [START_REF] Boffi | Finite element approach to immersed boundary method with different fluid and solid densities[END_REF] pour la méthode de la frontière immergée, et dans [START_REF] Burman | An unfitted Nitsche method for incompressible fluid-structure interaction using overlapping meshes[END_REF] pour une méthode de type Nitsche. L'incovénient majeur est, qu'avec ces approches, la stabilité ou la précision imposent des fortes restrictions sur le pas de temps (CFL parabolique) ou des itérations de correction. Dans cette note, nous présentons deux nouvelles méthodes numériques (Algorithmes 1 and 2) qui contournent ces difficultés. Leur caractère semi-implicite ou explicite dépend de l'ordre dans lequel les discrétisations spatiales et temporelles sont effectuées. Ces schémas généralisent (pour la première fois) les idées de [START_REF] Fernández | Incremental displacement-correction schemes for incompressible fluid-structure interaction: stability and convergence analysis[END_REF][START_REF] Fernández | Explicit Robin-Neumann schemes for the coupling of incompressible fluids with thin-walled structures[END_REF] au cas des maillages non compatibles. Les résultats des Propositions 2.1 et 3.2 et de la Section 4 montrent que la stabilité et la précision optimale ne sont pas soumises à des conditions CFL restrictives ou à des itérations de correction.

Introduction

We consider a fluid-structure interaction system in which the fluid is described by the Stokes equations and the structure is assumed to behave as a linear thin membrane or shell. The fluid domain is denoted by Ω f ⊂ R d (d = 2, 3). Its boundary is partitioned as ∂Ω f = Γ f ∪ Σ, where Σ represents both the fluid-structure interface and the solid domain (Figure 1, left). The exterior unit-vector normal to ∂Ω f is denoted by n. The coupled problem reads as follows: find the fluid velocity and pressure u :

Ω f ×R + → R d , p : Ω f × R + → R, the solid displacement and velocity d : Σ × R + → R d , ḋ : Σ × R + → R d such that            ρ f ∂ t u -divσ(u, p) = 0 in Ω f , divu = 0 in Ω f , u = 0 on Γ f , u = ḋ on Σ, (1)
     ρ s ∂ t ḋ + Ld = -σ(u, p)n in Σ, ḋ = ∂ t d in Σ, d = 0 on ∂Σ, (2)
complemented with standard initial conditions. Here, the constants ρ f and ρ s stand for the fluid and solid densities, respectively, while denotes the solid thickness. The fluid Cauchy-stress tensor is given by σ(u, p) def = -pI + 2µε(u) with ε(u) def = 1 2 ∇u + ∇u T , with µ denoting the fluid dynamic viscosity. At last, the abstract surface differential operator L describes the solid elastic effects.

Over the last years, significant progress has been achieved in the development and analysis of stable and accurate explicit coupling schemes (i.e., which invoke the fluid and solid solvers only once per time-step) for the approximation of the coupled problem (1)- [START_REF] Bukač | Fluid-structure interaction in blood flow capturing non-zero longitudinal structure displacement[END_REF]. All these studies (see, e.g., [START_REF] Burman | Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility[END_REF][START_REF] Guidoboni | Stable loosely-coupled-type algorithm for fluid-structure interaction in blood flow[END_REF][START_REF] Fernández | Incremental displacement-correction schemes for incompressible fluid-structure interaction: stability and convergence analysis[END_REF][START_REF] Fernández | Explicit Robin-Neumann schemes for the coupling of incompressible fluids with thin-walled structures[END_REF][START_REF] Burman | Explicit strategies for incompressible fluid-structure interaction problems: Nitsche type mortaring versus Robin-Robin coupling[END_REF][START_REF] Bukač | Fluid-structure interaction in blood flow capturing non-zero longitudinal structure displacement[END_REF][START_REF] Lukáčová | Kinematic splitting algorithm for fluid-structure interaction in hemodynamics[END_REF]) consider fitted meshes (Figure 1, left). For much applications, it is well known however that this assumption rapidly becomes cumbersome (e.g., large interface deflections, contacting structures, etc.). Within the unfitted mesh framework (Figure 1, right) splitting schemes of explicit nature are reported and analyzed in [START_REF] Boffi | Finite element approach to immersed boundary method with different fluid and solid densities[END_REF] using the finite element immersed boundary method, and in [START_REF] Burman | An unfitted Nitsche method for incompressible fluid-structure interaction using overlapping meshes[END_REF] for an unfitted Nitsche method. A major drawback of these approaches is that either stability or accuracy demands severe time-step restrictions (e.g., parabolic-CFL) or correction iterations. In this note, we present two new numerical methods (semi-implicit and explicit) which bypass these stability and accuracy issues. Their semi-implicit or explicit nature depends on the order in which the spatial and time discretizations are performed. These methods generalize (for the first time) the Robin-Neumann splitting paradigm of [START_REF] Fernández | Incremental displacement-correction schemes for incompressible fluid-structure interaction: stability and convergence analysis[END_REF][START_REF] Fernández | Explicit Robin-Neumann schemes for the coupling of incompressible fluids with thin-walled structures[END_REF] to the unfitted mesh framework. A priori energy and error estimates are stated for some of the variants. Their performance is illustrated via numerical experiments in a benchmark.

First discretize in space and then in time

The first class of methods introduced and analyzed in this note combines the unfitted mesh spatial approximation of [START_REF] Burman | An unfitted Nitsche method for incompressible fluid-structure interaction using overlapping meshes[END_REF] with the fractional-step time-marching of [START_REF] Fernández | Incremental displacement-correction schemes for incompressible fluid-structure interaction: stability and convergence analysis[END_REF][START_REF] Fernández | Explicit Robin-Neumann schemes for the coupling of incompressible fluids with thin-walled structures[END_REF].

Unfitted mesh spatial semi-discretization

Let {T s h } 0<h≤1 be a quasi-uniform family of triangulations of Σ. We consider a quasi-uniform family of meshes {T f h } 0<h≤1 which cover the fluid domain Ω f . We assume that each T f h is fitted to the boundary Γ f but not to Σ (Figure 1 tions associated to T s h and T f h will be denoted by X s h and X f h , respectively. For the approximation of the fluid we will consider the spaces

V f h def = v f h ∈ [X f h ] d v f h | Γ f = 0 and Q h def = X f h , while for the solid we set V s h def = v s h ∈ [X s h ] d v s h | ∂Σ = 0 .
The following notation will also be used

a f (u, p), (v f , q) def = 2µ ε(u), ε(v f ) Ω f -(p, div v f ) Ω f + (q, div u) Ω f and a s d, v s def = Ld, v s Σ ,
where the symbol (•, •) ω denotes the scalar product in L 2 (ω). Problem (1)-( 2) is approximated in space as follows: for t > 0, find

u h (t), p h (t), ḋh (t), d h (t) ∈ V f h × Q h × V s h × V s h , such that ḋh = ∂ t d h and    ρ f ∂ t u h , v f h Ω f + a f (u h , p h ), (v f h , q h ) + ρ s ∂ t ḋh , v s h Σ + a s (d h , v s h ) + s h (u h , p h ), (v f h , q h ) -σ(u h , p h )n, v f h -v s h Σ -u h -ḋh , σ(v f h , -q h )n Σ + γµ h u h -ḋh , v f h -v s h Σ = 0 (3) for all (v f h , q h , v s h ) ∈ V f h ×Q h ×V s h .
Here, γ > 0 denotes the Nitsche's penalty parameter and s h the stabilization bilinear form

s h (u h , p h ), (v f h , q h ) def = K∈T f h γph 2 µ ∇p h , ∇q h K + γ g µh ∇u h ∂K , ∇v f h ∂K ∂K
, where • ∂K denotes the jump across the boundary of the element K and γ p , γ g > 0 are free parameters. As shown in [START_REF] Burman | An unfitted Nitsche method for incompressible fluid-structure interaction using overlapping meshes[END_REF], the unfitted space semi-discrete formulation (3) is stable and delivers optimal first-order accuracy in the energy-norm.

Fully discrete formulation: semi-implicit coupling scheme with unfitted meshes

In what follows, τ > 0 denotes the time-step length, t n def = nτ for n ∈ N, ∂ τ x n def = 1 τ x n -x n-1 stands for the first-order backward difference and x denotes the r-th order extrapolation, that is, x = 0 if r = 0, x = x n-1 if r = 1 and x = 2x n-1 -x n-2 if r = 2. The spatial semi-discrete problem (3) is discretized in time with the incremental displacement-correction scheme displayed in Algorithm 1. Note that step [START_REF] Burman | An unfitted Nitsche method for incompressible fluid-structure interaction using overlapping meshes[END_REF] introduces and additional unknown, the intermediate solid velocity ḋn-

1 2 h
, which is implicitly coupled to the fluid through the solid inertial term. The remaining solid elastic contributions are treated explicitly (or ignored) in (4) via extrapolation. The end-of-step solid velocity ḋn h is obtained by solving the solid correction step [START_REF] Burman | Explicit strategies for incompressible fluid-structure interaction problems: Nitsche type mortaring versus Robin-Robin coupling[END_REF]. The implicit treatment of the solid inertia in (4) guarantees (added-mass free) stability, while the extrapolation of the solid elastic terms introduces a certain degree of fluid-solid splitting. Note that the intermediate solid velocity ḋn-1 2 h cannot be eliminated in (4) and, hence, the coupling scheme is not explicit. This is a major difference with respect to the case of fitted meshes considered in [START_REF] Fernández | Incremental displacement-correction schemes for incompressible fluid-structure interaction: stability and convergence analysis[END_REF][START_REF] Fernández | Explicit Robin-Neumann schemes for the coupling of incompressible fluids with thin-walled structures[END_REF].

Algorithm 1 Semi-implicit coupling scheme (i) Fluid with solid inertia sub-step: find

u n h , p n h , ḋn-1 2 h ∈ V f h × Q h × V s h such that      ρ f ∂ τ u n h , v f h Ω f + a f (u n h , p n h ), (v f h , q h ) + ρ s τ ḋn-1 2 h - ḋn-1 h , v s h Σ + a s (d h , v s h ) + s h (u n h , p n h ), (v f h , q h ) -σ(u n h , p n h )n, v f h -v s h Σ -u n h - ḋn-1 2 h , σ(v f h , -q h )n Σ + γµ h u n h - ḋn-1 2 h , v f h -v s h Σ = 0 (4) for all (v f h , q h , v s h ) ∈ V f h × Q h × V s h . (ii) Solid sub-step: find ḋn h , d n h ∈ V s h × V s h such that ḋn h = ∂ τ d n h and ρ s τ ḋn h - ḋn-1 2 h , v s h Σ + a s (d n h -d h , v s h ) = 0 (5) for all v s h ∈ V s h .
In the succeeding text, the quantities

E n h def = ρ f 2 u n h 2 0,Ω f + ρ s 2 ḋn h 2 0,Σ + 1 2 a s d n h , d n h and E n h def = ρ f 2 u n h -u(t n ) 2 0,Ω f + ρ s 2 ḋn h -ḋ(t n ) 2 0,Σ + 1 2 a s d n h -d(t n ), d n h -d(t n ) 1 2
respectively denote the discrete energy of the system and error at time t n . The symbol will indicates an inequality up to a multiplicative constant (independent of h and τ ). The following result, whose proof can be found in [START_REF] Fernández | Unfitted formulations and splitting schemes for fluid/thin-walled structure interaction[END_REF], states the stability and convergence properties of Algorithm 1. Proposition 2.1 (Stability and convergence) Let {(u n h , p n h , ḋn h , d n h )} n≥1 be given by Algorithm 1. For γ > 0 sufficiently large, there holds E n h E 0 h unconditionally for r ∈ {0, 1} and under the condition τ = O(h 6 5 ) for r = 2. Moreover, under the same conditions, there holds

E n h h + τ + τ 2 r-1 .
This result shows that Algorithm 1 preserves the stability and accuracy properties of the explicit coupling schemes introduced in [START_REF] Fernández | Incremental displacement-correction schemes for incompressible fluid-structure interaction: stability and convergence analysis[END_REF][START_REF] Fernández | Explicit Robin-Neumann schemes for the coupling of incompressible fluids with thin-walled structures[END_REF] with fitted meshes. Moreover, it overcomes the severe stability restrictions observed in [START_REF] Boffi | Finite element approach to immersed boundary method with different fluid and solid densities[END_REF] for the traditional time-marching schemes of the immersed boundary method. Algorithm 1 with r = 1 delivers unconditional optimal overall first-order accuracy. This is also significant progress with respect to the splitting schemes reported in [START_REF] Burman | An unfitted Nitsche method for incompressible fluid-structure interaction using overlapping meshes[END_REF], whose accuracy is known to be non-uniform in h.

First discretize in time and then in space

The main drawback of Algorithm 1 is its semi-implicit nature. In other words, step (4) has a computational complexity larger than a single fluid problem, due to the additional unknown ḋn-1 2 h

. In this section, we introduce a new explicit coupling scheme which overcomes this issue without compromising stability and accuracy. The fundamental idea consists in performing the space and time discretization reversely.

Robin-Neumann explicit coupling schemes

The explicit coupling schemes introduced in [START_REF] Fernández | Incremental displacement-correction schemes for incompressible fluid-structure interaction: stability and convergence analysis[END_REF][START_REF] Fernández | Explicit Robin-Neumann schemes for the coupling of incompressible fluids with thin-walled structures[END_REF] with fitted meshes, can be derived by applying the time splitting of Section 2.2 to the continuous problem (1)-( 2). This yields:

(i) Fluid with solid inertia sub-step: find

u n : Ω f ×R + → R d , p n : Ω f ×R + → R and ḋn-1 2 : Σ×R + → R d such that                    ρ f ∂ τ u n -divσ(u n , p n ) = 0 in Ω f , divu n = 0 in Ω f , u n = 0 on Γ f , u n = ḋn-1 2 on Σ, ρ s τ ḋn-1 2 - ḋn-1 = -Ld -σ(u n , p n )n on Σ. (6)
(ii) Solid sub-step: find

d n : Σ × R + → R d and ḋn : Σ × R + → R d such that ḋn = ∂ τ d n and    ρ s τ ( ḋn - ḋn-1 2 ) + L(d n -d ) = 0 in Σ, d n = 0 on ∂Σ. (7)
Contrarily to Algorithm 1, the intermediate solid velocity ḋn-1 2 in the fluid step (i) can be eliminated through [START_REF] Caiazzo | Projection schemes for fluid flows through a porous interface[END_REF] 4 . This yields the following explicit Robin-Neumann time splitting on the interface:

σ(u n , p n )n + κu n = κ ḋn-1 + g on Σ, ρ s ∂ τ ḋn + Ld n = -σ(u n , p n )n on Σ, (8)
with

g def = ρ s ∂ τ d + σ(u , p )n and κ def = ρ s /τ .

Fully discrete formulation: explicit coupling scheme with unfitted meshes

The fundamental idea consists in performing directly an unfitted interface treatment (à la Nitsche) of the Robin-Neumann time splitting [START_REF] Fernández | Unfitted formulations and splitting schemes for fluid/thin-walled structure interaction[END_REF], by extending the arguments introduced in [START_REF] Burman | An unfitted Nitsche method for incompressible fluid-structure interaction using overlapping meshes[END_REF][START_REF] Caiazzo | Projection schemes for fluid flows through a porous interface[END_REF]. The proposed fully discrete scheme is based on the following result (see [START_REF] Fernández | Unfitted formulations and splitting schemes for fluid/thin-walled structure interaction[END_REF] for a proof). Proposition 3.1 (Consistency) Let {(u n , p n , ḋn , d n )} n≥1 be given by ( 6)- [START_REF] Fernández | Incremental displacement-correction schemes for incompressible fluid-structure interaction: stability and convergence analysis[END_REF]. Then, there holds [START_REF] Fernández | Explicit Robin-Neumann schemes for the coupling of incompressible fluids with thin-walled structures[END_REF] is the fact that for κ → ∞ (i.e. τ → 0) we formally retrieve the unfitted formulation [START_REF] Burman | Stabilization of explicit coupling in fluid-structure interaction involving fluid incompressibility[END_REF]. Taking successively v s h = 0 and (v f h , q h ) = (0, 0) in ( 9) we obtain the fully discrete method reported in Algorithm 2. Note that the resulting coupling scheme is explicit. The following result, whose proof can be found in [START_REF] Fernández | Unfitted formulations and splitting schemes for fluid/thin-walled structure interaction[END_REF], guarantees the energy stability of Algorithm 2 for r = 0. Proposition 3.2 (Stability) Let {(u n h , p n h , ḋn h , d n h )} n≥1 be given by Algorithm 2 with r = 0. For γ > 0 sufficiently large, we have

             ρ f ∂ τ u n , v f h Ω f + a f (u n , p n ), (v f h , q h ) + ρ s ∂ τ ḋn , v s h Σ + a s (d n , v s h ) + γκµ γµ + κh u n - ḋn-1 , v f h -v s h Σ - κh γµ + κh σ(u n , p n )n, v f h -v s h Σ + u n - ḋn-1 , σ(v f h , -q h )n Σ - h γµ + κh σ(u n , p n )n, σ(v f h , -q h )n Σ - γµ γµ + κh g , v f h -v s h Σ + h γµ + κh g , σ(v f h , -q h )n Σ = 0 (9) for all (v f h , q h , v s h ) ∈ V f h × Q h × V s h . A salient feature of
E n h ≤ E 0 h . Algorithm 2 Explicit coupling scheme (i) Fluid sub-step: find u n h , p n h ∈ V f h × Q h such that                ρ f ∂ τ u n h , v f h Ω f + a f (u n h , p n h ), (v f h , q h ) + s h (u h , p h ), (v f h , q h ) + γκµ γµ + κh u n h - ḋn-1 h , v f h Σ - κh γµ + κh σ(u n h , p n h )n, v f h Σ + u n h - ḋn-1 h , σ(v f h , -q h )n Σ - γµ γµ + κh g h , v f h Σ - h γµ + κh σ(u n h , p n h )n, σ(v f h , -q h )n Σ + h γµ + κh g h , σ(v f h , -q h )n Σ = 0 for all (v f h , q h ) ∈ V f h × Q h . (ii) Solid sub-step: find ḋn h , d n h ∈ V s h × V s h such that ḋn h = ∂ τ d n h and      ρ s ∂ τ ḋn h , v s h Σ + a s (d n h , v s h ) + κh γµ + κh σ(u n h , p n h )n, v s h Σ - γκµ γµ + κh u n h - ḋn-1 h , v s h Σ + γµ γµ + κh g h , v s h Σ = 0 for all v s h ∈ V s h .
The stability of Algorithm 2 with r = 1, 2 as well as the convergence properties of all the variants are illustrated in Section 4 via numerical experiments.

Numerical experiments

In order to highlight the stability and accuracy of the proposed schemes, we consider the numerical example of [START_REF] Burman | An unfitted Nitsche method for incompressible fluid-structure interaction using overlapping meshes[END_REF]. We compare the results obtained with Algorithms 1 and 2 with those obtained with a first-order fully implicit scheme using fitted meshes. The convergence histories are displayed in Figure 2. In order to highlight the uniformity of the convergence in h, we have refined both in time and in space at the same rate, τ = O(h). In spite of their different semi-implicit and explicit nature, Algorithms 1 and 2 deliver practically the same behavior: stability is obtained with all the variants, optimal first-order convergence is obtained with the extrapolated variants (r = 1, 2, unfitted) and the implicit (fitted) scheme while a sub-optimal convergence rate is exhibited by the non-extrapolated ones (r = 0, unfitted).

Figure 1 .

 1 Figure 1. (left) and unfitted (right) fluid meshes.

Figure 2 .

 2 Figure 2. Time convergence history of the solid displacement in the relative elastic energy norm using Algorithm 1 (left) and Algorithm 2 (right) with τ = O(h).

(Mikel Landajuela).

1 . Work funded by the french National Research Agency (ANR) through the EXIFSI project (ANR-12-JS01-0004).


Dominique Vaufreydaz
email: dominique.vaufreydaz@inria.fr

Amaury Nègre
email: amaury.negre@imag.fr

MobileRGBD, An Open Benchmark Corpus for mobile RGB-D Related Algorithms

Since the commercialization of low cost RGB-D sensors, like the Kinect, more and more indoor robots have been equipped with this kind of sensors to perform tasks as people tracking or gesture recognition. Nevertheless, as far as we know from the literature, studies do not consider the limits of the sensors in term of motion speed, position of the sensor on the robot, etc. In this work, we propose to provide a corpus dedicated to low level RGB-D algorithms benchmarking. Originality of our approach is the use of dummies in order to play static users in the environment. This idea let us vary other variables that can impact algorithm performance: linear/angular speed of the robot, trajectory of the robot, RGB-D sensor height and vertical angle of view, number and relative position of dummies and furniture position. This paper first describes the experimental platform used to perform the acquisitions and the environment setup required to reproduce the dataset. Then, a precise description of all available data is given. We will see that, as this corpus contains a lot of configurations, it will allow researchers to investigate how these variables impact the results of their algorithms.

I. INTRODUCTION

Body shape, gesture, body language detection and evaluation are not new areas of research. Several years ago, one can use computer vision algorithms in order to detect pedestrians, bodies or objects from a video stream. Later on, researchers have used Time-Of-Flight (TOF) cameras in order to achieve this task but these devices are expensive. Nowadays, many researches on this topic are done using Kinect of other RGB-D sensors [START_REF] Cruz | Kinect and rgbd images: Challenges and applications[END_REF].

Many robotic researches, especially but not only in researches about companion robot, are using RGB-D sensors like the Kinect for many tasks: 3D SLAM, people tracking and skeleton tracking, gesture recognition, body pose estimation, engagement toward a robot or for elderly at home, detecting falls or searching for a fallen person ( [START_REF] Sturm | A benchmark for the evaluation of rgb-d slam systems[END_REF], [START_REF] Kim | Gesture recognition for moving rgb-d sensor[END_REF], [START_REF] Mastorakis | Fall detection system using kinect's infrared sensor[END_REF], [START_REF] Negin | A decision forest based feature selection framework for action recognition from rgb-depth cameras[END_REF], [START_REF] Buys | An adaptable system for rgb-d based human body detection and pose estimation[END_REF], [START_REF] Benkaouar | Multi-Sensors Engagement Detection with a Robot Companion in a Home Environment[END_REF], [START_REF] Anzalone | Multimodal people engagement with icub[END_REF]). We totally agree with [START_REF] Berger | The role of rgb-d benchmark datasets: an overview[END_REF], benchmarking is one pillar of research in these domains: how can one states that an algorithm performs better without a common reference? For instance, in [START_REF] Davari | Automatic fall detection for elderly by using features extracted from skeletal data[END_REF], a fall detection system is depicted with an interesting approach and good performances but description of the train/test corpus is not detailed enough. Thus, this corpus is not reproducible and results not comparable to other approaches.

RGB-D related algorithms are difficult to benchmark as there is no reference corpora with ground truth labels for every task or/and condition. Focusing on all body related tasks, it is even more difficult to record data. Labeling depth data with skeletons is not an easy task and recording many times people enhanced variability but increase difficulty of the labeling task.

Several corpora using RGB-D are available for many tasks. One can find corpora for language signing [START_REF] Stefanov | A kinect corpus of swedish sign language signs[END_REF] or mental state computation [START_REF] Mahmoud | 3d corpus of spontaneous complex mental states[END_REF] for instance. [START_REF] Wolf | The liris human activities dataset and the icpr 2012 human activities recognition and localization competition[END_REF] proposed a Human activities dataset for ICPR 2012 human activities recognition and localization competition. This dataset exposes humans in everyday-life tasks recording using the Kinect sensors. Labeling is done on human activities, i.e. at a high level of annotation. One cannot state performances of skeleton or face detection using this corpus even these features can be used for the activity recognition. Few corpora are available on with a mobile robot equipped with an RGB-D sensor. In [START_REF] Sturm | A benchmark for the evaluation of rgb-d slam systems[END_REF], the RGB-D data are dedicated to SLAM. Mobile paths were done using a man handled and a robot mounted Kinect device. For the robot recordings, linear speed is from 0.1 to 0.23 m.s -1 and angular speed is around 12 deg.s -1 . Ground truth labels are provided by a calibrated motion capture system, but no body features are provided. One can find a more correlated approach in [START_REF] Kim | Gesture recognition for moving rgb-d sensor[END_REF]. The authors use a Kinect mounted on a cleaning robot in order to evaluate upper body skeleton tracking. People were asked to play predefined gestures in front of the robot while it is rotating. In this study, we have no information about the robot or its angular speed. Moreover, variability in this case is not only due to the speed of the robot but also to human gestures. Between two records, humans' moves are different, thus can introduce variability on detection algorithms. These artifacts can lead to different results for several speeds, not depending on the speed itself but on external variables.

In this paper, we first present our design path and objectives in section II. Section III depicts our experimental environment, i.e. our experimentation room and robotic platform. We give clues about control of the robot in order to make reproducible recording trajectories and RGB-D sensor vertical angle calibration. Gathered data using the Kinect 2 device mounted and the robot sensors are detailed in section IV.

II. DESIGN PATH AND OBJECTIVES

From our experience of using RGB-D sensors on mobile robot, we know it could be difficult to isolate which changes/variables impact performance of a specific algorithm: is it due to the robot speed, the number of persons at the same time, their relative position, etc. Sometimes, we even improve performance by changing position or field of view of the mounted RGB-D sensor on the robot. Looking at design of some mobile platforms like the Hobbit [START_REF] Fischinger | Hobbit-the mutual care robot[END_REF], the turtlebot 1 [START_REF] Xiong | An exclusive human-robot interaction method on the turtlebot platform[END_REF], a tobotic Wheelchair [START_REF] Vasquez | Human aware navigation for assistive robotics[END_REF] or our version of the Kompaï [START_REF] Benkaouar | Multi-Sensors Engagement Detection with a Robot Companion in a Home Environment[END_REF], positions of using mounted RGB-D sensor differ. Most of the time, it was parallel to the floor plane. These choices are often driven by human, technical or design considerations without evaluating impact on the algorithms. From all the reasons above, we can extract a set of variables that we want to address in our corpus:

• linear and angular robot speed;

• number and position of persons in the environment;

• furniture placement and body occlusions;

• height and vertical field of view of the RGB-D sensor, i.e. looking down or up at several positions;

• body orientation toward the sensor;

We thought to reverse the corpus recording paradigm. Our goal is to facilitate ground truth annotation and reproducibility of records among speed, trajectory and environmental variations. As we want to get rid of unpredictable human moves, in our benchmark corpus, we use dummies (see 2). Interest of dummies resides in the fact that they do not move between two recordings. It is possible to record the same robot move in order to evaluate performance of detection algorithms varying speed.

Our goal is to provide a benchmark corpus for "low level" RGB-D algorithm family like 3D-SLAM, body/skeleton tracking or face tracking using a mobile robot. Using this open corpus, researchers can find a way to answer several questions:

• what is the algorithm performance in multiples conditions?

• on a mobile robot, what is the maximum linear/angular speed supported by the algorithm?

• which variables impact the algorithm?

• evaluate suitable height/angle of the mounted RGB-D sensor to reach goals: monitoring everyday live is different from searching fallen persons on the floor;

• finally, what is the performance on an algorithm with regards to others?

III. EXPERIMENTAL SETUP

A. Room description

The experimental room is a flexible space designed to be representative of a home-like environment. It has an 'L' shape with a kitchen place with a sink, a diner and a lounge space. The size of the room is 6x8 m. A view of the kitchen space is shown on the fig. 2. Except the sink, furniture can be moved ad libitum to reflect experimental needs.

In our setup, dummies are in the kitchen. The robot moves from the dining room toward the dummies. Some furniture is available to create different setups and placed to make some occlusions: sofa, table and chairs. Augmented reality tags are on the walls and let us have more 3D information that information gathered from the robot laser range finder (see III-B).

B. Platform description

Our recording platform is depicted on fig 1. For the robotic part, we used a Robulab10 2 mobile platform with 2 lateral propulsive and 2 castor wheels. Each propulsive wheel can been control independently. The Robulab10 is equipped with a laser range finder, 8 ultrasound and 16 infrared telemeters. On its top, we mounted a flexible structure with a Kinect 2. This structure lets us change manually the height of the sensor and, using servomotor, one can programmatically modify the vertical angle of view. Last, a laptop is used for control and record facilities.

C. Localization and trajectory following

The robot trajectory repeatability is one of the key aspects to compare different configurations. In this sense, an absolute accurate robot localization is required. To provide such localization a map of the room is exploited in combination of the laser sensor data. The chosen solution consists on a variant of Iterative Closest Point algorithm [START_REF] Censi | An icp variant using a point-to-line metric[END_REF] based on a point-to-line measure to compute the transformation that best matches the laser point cloud to the walls and static furniture (see Fig. 3). This localization is performed in real time and we experimentally measured a centimeter accuracy. Once the robot localization is known, we implement a rectilinear trajectory following algorithm. The robot command is decomposed on a linear speed and an angular speed. To control the angular speed, we first compute a control point in front of the robot, located in the linear trajectory (see Fig. 4). Then we calculate the error (ψ) as the angle between the robot direction and the control point direction. This error is injected to a PID controller to obtain a filtered angular speed command. For the linear speed control, we simply use a trapezoidal speed profile, decomposed in three steps to avoid hard acceleration. In the first step, the speed v increased linearly with the time t:

v(t) = a • t 2
In the second step, the speed is constant at requested value. Last, the robot decelerates progressively, the speed decrease is computed using the square root of the distance to the goal d:

v = √ 2 • a • d

D. Sensor height and tilt mechanism

Height of the RGB-D device, i.e. its basement position relatively to the floor, can be modified manually. This position varies from 48 to 128 cm. The sensors are 4 cm over the basement, thus, the recording positions of the mounted Kinect 2 are from 52 to 132 cm with a 10 cm step. In our hardware structure, we added a servomotor to rotate vertically the Kinect. Behind this idea, we want to investigate if it is more suitable to have a lower Kinect looking up, a central or a higher one looking straight or down on a future robot. This setup also deserves to improve search algorithm for a fallen person, i.e. we must have a floor view using the Kinect. From the height of the RGB-D sensor (h), computation of the percentage of a 1.8m person in the field of view at l=4m for a specific angle ka. Z1 and Z2 are respectively lowest and highest visible points at a 4m distance. One can find an example in table I.

Step of the servomotor is (360/4096) • < 0.1 • . To calibrate it for our specific angle set, we first extracted from the depth data the normalized coefficients of the floor plane x, y, z and w the actual height of the sensors. For each servomotor position, we computed the Kinect angle (ka) in regards to the floor: ka = atan2(F loorP lane.z, F loorP lane.y)

Obviously, recording angles depend on the height of the RGB-D sensor. On lower position, we did not record floor view and on higher position we did not record ceil view. As the max depth value is 4.5m, We decided to record every angle that permits at 4m to see at least 50% of a 1.80m body (see Fig 5). Using ka, we computed Z1 and Z2, respectively the lowest and highest visible points at a 4m distance. The percentage of visible body is then computed directly: visible % = ||Z1, Z2 ∩ 0, 1.80m||/1.80 Fig. 6. Varying the height of the sensor (h in fig. 5), compute percentage of a 1.80m body in the field of view at l=4m regards to a set of angles. For each height, retained recording angle are over 50%. Fig 6, shows all body percentage values varying h. In the table I, one can find example using the lowest sensor height (h=0.52m). In this case, we recorded using angles from -20 • to 35 • using a 5 • step. For the highest position (h=1.32m), retained angles are from from -35 • to 20 • with the same step.

IV. CORPUS ACQUISITION

A. Scenarios

Among variables, our scenarios first must handle dummies and furniture placements within the experimental space. Dummies are set by one or two in several positions: in front of the rear wall, in front of the sink, sited on the sofa or behind a table to create partial body occlusion for instance (see fig. 2). Other variables for our recording scenarios are recording speed and path. On fig. 7, one can see some example of recording trajectories (more trajectories are tackled in the corpus). Each blue arrow represents a robot trajectory recorded at several speeds forward and backward. We do strait trajectories ending far or closer to the rear wall. Curved paths and pure rotations are also used. Last, as seen on fig. 7, trajectory at 45 • are followed by our robot.

For each trajectory, several records are made setting the linear speed from 0.1 to 1.1 m.s -1 with a 0.1 step. The maximum linear speed of the propulsive wheels 1.3 m.s -1 , thus we must limit the maximal speed in order to be able to make the robot turn. For the pure rotations, angular speed is set from 0.1 rad.s -1 to 2 rad.s -1 . We envision to record more complex paths among a full home-like setup and in some corridors.

B. Data acquisition

This project started with the release of the new Kinect 2 sensor 3 . At the recording time, OpenNI under Linux 4 does not handle correctly this new device. Thus, we choose to use the sensor under the Windows Kinect SDK. We recorded synchronously all available streams from our mobile platform (see section III-B). Features are all robot centered.

The Kinect sensor exposes several data streams. Some of them are shown in fig. 1. Some difficulties raised when realtime acquisition of all the sensors equipping the Kinect is done. Asking for multiple streams at the same time can lead to variations on streams frame rate. For instance, when a skeleton is detected, the RGB frame rate slightly decreases. Optimization on the acquisition code was made in order to record all streams without information loss. An equal care was made to gather data from the Robulab10 platform.

Finally, we record at the highest possible frame rate several robot centered information. RGB, depth and infrared streams, detected skeletons are gathered from the Kinect 2. From the robot, laser range finder, acoustic and infrared telemeters, battery level, wheel odometry, current linear and angular speeds are monitored. All information is stored in uncompressed format in order to remove coding/decoding artifacts on video streams for instance. Each frame/event is tagged using timestamps in ms from epoch time 5 and can be synchronized for viewing or processing. The data recorded to build the corpus are summarized in the table II.

1) Videos streams:

The resolution of the RGB image is of 1920 × 1080 pixels at 30 frames per second. The RGB sensor has a 70 • horizontal and 60 • vertical field of view wide-angle lens. Data are stored in 2 files. The first one is a binary raw concatenation of all YUV received images. The second one is a text file containing all timestamps associated with video frames. A compressed video file is provided for human needs. Portable C++ source code is also provided for reading frames using OpenCV6 .

2) Depth and Infrared streams: In the Kinect 2, the depth stream is a real Time-Of-Light (TOF) approach whereas previous Kinect versions that use structured light. It is also an active IR device, thus there is an IR video stream associated to the depth information. The depth range of the Kinect is limited from 0.4 meter to 4.5 meters. The depth is measured in meters from the camera along the Z axis, X and Y axis are in pixel coordinates. The resolution of the depth and IR images is of 512 × 424 pixels. The frame rate is at maximum 30 frames per second. As for the RGB streams, IR and depth streams are both stored using a binary raw file and associated timestamps in a separate text file. 3) Body information: The Kinect on Windows SDK supports up to 6 skeletons tracked at the same time (it used to be 2 with the previous sensor version). Each skeleton is associated to a bag of pixels from the depth image. They now contains 25 joints (20 previously) more morphologically placed within the body shape. New joints concern neck, thumbs and hand tips. 3D information about joints is enhanced by quaternions giving rotations of body parts. We also have information about face detected in the RGB stream and/or from the depth stream. For portability, timestamped body data are saved in JSON format 7 in a text file. A portable C++ source code is provided for reading data.

4) Robulab10 sensors: As described, Robulab10 is equipped with many sensors. Even if we use the laser data for controlling the robot (see III-C), the nine ultrasound and the sixteen infrared telemeters are monitored. Wheel odometry, internal Robulab10 localization, commands sent to the robot and battery level are also gathered. All these data are saved in JSON with timestamps.

C. Labels

To make reproducible experimentation, one needs ground truth information and labels. For each experimental setup, several information are provided:

• a static map and a SLAM Based map of room and furniture;

• a 2D position and a 3D bounded box for each dummy;

• 10 seconds of full recording of each dummy are available (see section IV-B). Among this information, one can find an automatic skeleton from the Kinect and its projected information in the room space. This recording is done with a localized and still robot, using a Kinect at middle height and parallel to the ground. The robot is manually positioned at a strategic place, 7 http://www.json.org/ i.e. a place where we have an skeleton for the dummy and at least 2 visual tags (see Fig. 2).

• Another skeleton is manually annotated using the Depth data.

D. Corpus availability

The corpus is freely available for research teams through a web site 8 . As we said, it will come with C++ source code provided to read synchronously data. These source codes work under Window and Linux. We envision providing also a ROS "bag" 9 in order to facilitate integration of our corpus. After the final release of the Kinect 2, we plan to release our recording source code.

We have several hundreds of Gigabytes available for the community. We are preparing a web solution to conveniently distribute it. Indeed, we need to find a solution to distribute it the right way. We plan to propose a website where people can download a selection of the available data. Doing this, people may choose first subsets to achieve first experiments. Then, they can better choose other subsets in order to validate their first results.

We encourage researchers from other research teams to contribute to our effort. Using their own robot, they can use the same techniques to record other scenarios in other conditions.

V. CONCLUSION

In this paper, we have presented MobileRGBD a new freely available dataset for benchmarking RGB-D related algorithms on a mobile platform. At the writing time of this paper, we are still recording data. This corpus contains color images, depth maps and IR images, body information from the Kinect 2 and localization information of our Robulab10 platform. Originality of our corpus is the use of dummies in order to play static users in the environment. This idea let us vary other variables that can impact algorithm performance: linear/angular speed of the robot, trajectory of the robot, RGB-D sensor height and vertical angle of view, number and position of dummies and furniture position.

We propose a dataset that allows researchers to evaluate what is the performance of their algorithms and which variables impact their results. Knowing that, they can make more enlighten choices for the design of a future robotic platform or to solve specific problems. A multi-platform specific C++ source code is providing to facilitate the use of our data. The corpus will be available though a Web Site http://www-prima. inrialpes.fr/Vaufreydaz/OpenRGBDBenchmarkCorpus/.

  Fig. 1. Robotic platform used for recording: a Robulab10 mobile robot from the Robosoft compagny. It is equipped with a laser ranger finder, a Kinect 2.0 and a laptop computer.

Fig. 2 .

 2 Fig.2. View of the experimental room. The room is designed to simulate a real living apartment. Some augmented reality tags have been placed around the two dummies.

Fig. 3 .

 3 Fig.3. Robot localization, an ICP point-to-line method is used to compute the transformation that minimize the distance between laser points and map.

Fig. 4 .

 4 Fig. 4. Trajectory servoing.

Fig. 5 .

 5 Fig. 5.From the height of the RGB-D sensor (h), computation of the percentage of a 1.8m person in the field of view at l=4m for a specific angle ka. Z1 and Z2 are respectively lowest and highest visible points at a 4m distance. One can find an example in table I.

Fig. 7 .

 7 Fig. 7. Exemples of robot moves while recording. The dummies (green rectangles) and furniture (red lines) setup is the one presented on figure 2. Blue arrows represent robot trajectories.

Fig. 8 .

 8 Fig. 8. Robot sensor view. At left, the laser range finder view. In the middle, the RGB stream. Last, the infrared view. The depth view (not on this figure) is aligned with the infrared one.

TABLE I .

 I KNOWING THE HEIGHT OF THE SENSOR (HERE h=0.52m), COMPUTATION OF Z1, Z2 (SEE FIG. 5), AND PERCENTAGE OF BODY IN THE FIELD OF VIEW FOR A SET OF ANGLES. IN BOLD, RECORDING ANGLES WITH PERCENTAGE OVER 50%.

  with 25 joints and body parts rotations and faces Skeletons computed using OpenNi and face detection using OpenCV will be provided.

	Data		Sensor	Information	Frame rate
	Telemeters distances	Laser range finder	20 meters maximum	12.5Hz
	Ultrasound distances	Ultrasound telemeters	3 meters maximum	12.5Hz
	IR distances	IR telemeters	1.5 meters maximum	12.5Hz
	Commands	Control program	linear and angular speeds, stop command	12.5Hz
	Odometry	Robulab10	internal localization information	12.5Hz
	Battery level	Robulab10	in percentage	1Hz
	Body		Kinect 2	maximum 6 skeletons 30Hz max
	RGB Video	Kinect 2	1900x1080p	30Hz
	Depth Video	Kinect 2	512x424p from 0.4 to 4.5m	30Hz
	IR Video	Kinect 2	512x424p	30Hz
	TABLE II.	SUMMARY OF COLLECTED DATA IN THE CORPUS WITH THEIR CORRESPONDING SENSORS, INFORMATION AND MAXIMAL ACQUISITION
				FRAME RATE.

The Robulab10 platform is sold by the Robosoft compagny http://www. robosoft.com/

We are members of the Kinect 2 for Windows beta test program.

http://www.openni.org/?

The timestamps are expressed in time since 1 st of January, 1970.

http://http://opencv.org/

http://www-prima.inrialpes.fr/Vaufreydaz/OpenRGBDBenchmarkCorpus/

ACKNOWLEDGMENT

The authors would like to thank Inria and French Ministry of Education and Researches for their support. This work was done using the Amiqual4Home facilities (ANR-11-EQPX-0002).


Adolphe Kimbonguila

Mohamed Guessasma
email: mohamed.guessasma@u-picardie.fr

Patrice Coorevits

Jérôme Fortin

Emmanuel Bellenger

Numerical modeling of granular material fragmentation by using Discrete Element Method

Keywords: granular material, DEM, fragmentation, strength criteria

The aim of this study is the modeling of the fragmentation of granular material subjected to compressive loading by using Discrete Element Method (DEM). The proposed model depends mainly on two mechanical aspects which are: the mechanical behavior of the granular material and the yield strength predicting the fragmentation of grains. In order to better understand the mechanical phenomena, we have studied the contact force chains according to the granulometry of the material and the fragmentation modes of grains. The proposed modeling has been tested on a numerical simulation of oedometric test by using our software MULTICOR.

Introduction

The crushing (or grinding) process is usually used in many industrial fields such as: pharmaceutical, mineral or food. It is well known that this process remains economically expensive, essentially due to the energetic cost and the low quality of finished products (heterogeneity of grain size and shape). In fact, the economic stakes are large as for the optimization of such process. Two solutions are commonly proposed to cure these problems: the first solution is based on the technological aspects while the second one is based on physical aspects. As regards the technological solutions, they are mainly based on the choice of the engine type to be used in order to obtain a powerful installation of crushing process and which could be adapted for granular material. However, these technological optimizations haven't given satisfactory results because the production cost remains so much expensive. In fact, it seems to be important to take into account physical aspects in order to better understand the mechanical phenomena during the process. Few years ago, many authors have suggested predictive models about physical parameters (friction and wear coefficients). Generally, these models are based on approaches of chemical engineering, but they are not appropriated to study the force chains distribution into the granular material and the prediction of fragmentation process of grains. Our goal, through out this study, consists to develop a discrete model for granular material behavior subjected to compressive stresses by using DEM [START_REF] Aström | Fragmentation of grains in a two-dimensional packing[END_REF], (Morris et al. 1998) and [START_REF] Morikawa | Analytical studies on local damage to reinforced concrete structures under impact loading by discrete element method[END_REF]. This modeling depends on two mechanical aspects which are: the mechanical behavior of the granular material and the yield strength for the prediction of grain fragmentation. In order to better understand the mechanical phenomena, we must answer the following questions:

How the contact forces are distributed into the granular material subjected to compressive stresses and how the contact forces are organized according to the granulometry of the material?

Which mechanisms are brought into play in the fragmentation of grains and which impact the granulometry have on the fragmentation modes?

Which overall behavior for the granular media during the process of crushing and which contribution of each particle during the process of fragmentation?

Mechanical approach with DEM

The conventional DEM allows to model really deformable particles as well as complex shapes (from the ellipsoid to the polygon). Here, we have studied the simple case of nondeformable and non penetrable particles in 2D with the computational software MULTICOR, developed by Fortin et al. (2004). The coordinates and the rotations of Euler are the configuration parameters q. The gyroscopic and centrifugal terms are equal to zero in 2D. The matrix of generalized mass M of the system doesn't depend on q that is diagonal block. The mechanical equation ( 1) can be written in the following form:

( ) α R t , q , q F q M ext + = & & & (1)
where F ext represents the known external forces and α R the unknown interior forces related to contact reactions and α the contact number.

In a system composed of p heterogeneous particles (Figure 1), the critical parameter for the modeling time is the maximum number of interactions between particles. The more the interaction range is important the more we have to test the possible interactions between particles. MULTICOR uses the partitioning method coupled to a connectivity table (Fortin and Coorevits, 2003). This technique allows to reduce considerably the computational time. In that case, the computational time no longer increases like O(p 2 ) but only like O(p), which is almost optimal. For each couple of particles Ω i and Ω j which may enter in contact, is associated with a local reference whose axes are oriented according to the two unit vectors n and t, respectively normal and tangential vectors in the contact plan.

The normal n is directed from Ω i and Ω j . The variables put in duality are ij u & , the relative local velocity of Ω i compared to Ω j , and the contact reaction ij r of Ω i and Ω j . In the local base, they are written by equations (2): Modelizing grains as rigid solids in punctual contact for each grains p occupying the volume V p , the expression of mean stress tensor is given by the following equation ( 4):

n r r r , n u u u ij n ij t ij ij n ij t ij + = + =
( )

        ⊗ + - ⊗ = ∫ ∑ p V c p p r q dV q g q V 1 α α α ρ σ & & (4)
where g denote the gravity, q the position vector, q& & its acceleration and α r the contact reaction acting on the particle at the point of contact α c q . According to the additivity on internal moments in the sense of [START_REF] Moreau | Numerical investigation of shear zones in granular materials[END_REF], the mean stress tensor σ on the REV can be written as follows (5):

∑ ∈ = V p p p V V 1 σ σ (5)
where e denotes the void ratio.

In [START_REF] Fortin | Construction of an average stress tensor for a granular medium[END_REF] and [START_REF] Guessasma | Numerical modelling of mechanical tests using discrete element method[END_REF], the authors have shown the influence of dynamic effects on the mean stress tensor. Finally, the tensor p σ is composed in three parts: contact reaction stress, gravity stress and stress part due to the inertial forces. Accordingly, the contribution of inertial forces insures the symmetry of the stress tensor in the sense of Cauchy.

Modeling of fragmentation of granular material submitted to compressive stresses

The mean stress tensor presented in the previous section will be used in the modeling of fragmentation of granular material. The distribution of force chains into the granular material resulting from the compressive loading induces stress fields which are responsible of the fracture of grains. To predict the fragmentation of grain, we have considered the Mohr Theory of Failure, also known as the Coulomb-Mohr criterion or internal-friction theory (figure 3), where c σ and t σ are respectively the uniaxial tensile and compression strength. This theory is used in predicting the failure of brittle materials, and is applied to cases of plane stresses.

Figure.3. Coulomb-Mohr criterion

The Coulomb-Mohr criterion states that failure occurs when the two principal stresses are outside the criterion domain in the principal stress plane. Several experimental tests using granular materials, under compressive loading, have been carried in order to understand the different fragmentation modes (figure 4a). In our case, we have considered the fragmentation mode of grains which is represented in figure 4a [START_REF] Tsoungui | Partial pressure supported by granulometric classes in a polydisperse packing[END_REF]. The corresponding numerical model is composed of 12 grains resulting from the subdivision of the initial grain. Let us consider a i the radius of a given grain. We suppose that the Coulomb-Mohr criterion is not verified by the principal stresses into the grain for a given loading. Accordingly, the fragmentation mode of grain will be performed according to the following subdivision: the initial grain is divided in 4 gains with radius equal to a i /2 and a i /3, and 8 gains with radius equal to a i /6 and a i /11 (figure 4b).

Numerical simulation of oedometric test

In this part, we investigate a numerical two-dimensional model for the fragmentation of a granular material under pressure (Figure 5). The calculations are carried out with MULTICOR software. We consider that the behavior of granular material is perfectly rigid with density equal to 1000kg/m 3 and friction coefficient about 0.2. The granular packing contains 32 grains having diameter equal to 10mm and uniaxial compression strength about 5MPa. The granular system is subjected to a pressure velocity about 0.5m/s.

Conclusion

In conclusion, we have shown that the contact force chains are crucial to the fragmentation of granular material under pressure. By considering the gravity forces after fragmentation, the rearrangement and distribution of grains seems more In addition, the overall behavior of granular material reproduces phenomena which are observed in granular media like vault effect or dilatancy.

Figure. 2 .

 2 Figure.2. REV of granular material.

Figure. 4 .

 4 Figure.4. (a) Fragmentation modes and (b) Numerical model.

Figure. 5 .

 5 Figure.5. Oedometric test of granular packing.

Figure. 6 .

 6 Figure.6. (a) Forces chains distributions and (b) Uniaxial compression stress.

& & &

The introduction of Coulomb's friction µ leads to a non linear problem which can't be solved by a linear programming method. Unlike the usual approach, the bipotential method leads to a single variational principle and an inequality (Fortin et al. 2004). By using Usawa's algorithm, we obtain a resolution algorithm of the constitutive law based on the predictive-corrective scheme expressed by equations (3):

where γ is a numerical parameter and µ K the Coulomb's cone.

The projection of the prediction ij τ on the Coulomb's cone µ K leads to one of these states: non contact, contact with adherence or sliding contact. 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. The convergence is obtained when the force confirms the unilateral contact law with dry friction.

Mean stress tensor of granular material

In order to understand the different mechanisms governing the fragmentation of the grains subjected to compressive loading, we must be able to calculate the mean stress tensor into the granular material. The knowledge of stress field, acting in the granular material, permits to give a prediction of fragmentation of grains with respect the yield strength of the granular material. Let us consider a Representative Elementary Volume (REV) containing N p grains separated by voids (Figure 2).


Adolphe Kimbonguila

Mohamed Guessasma

Patrice Coorevits

Jérôme Fortin

Granular Crushing Of

Media

CRUSHING OF GRANULAR MEDIA BY A DISCRETE NUMERICAL MODELING

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The solutions commonly suggested in order to "optimize the performances" of crushing process are based on the technological or physical aspects. As regards the technological solutions, they are mainly based on the choice of the engine type to be used to obtain a powerful installation of crushing process. Unfortunately, these technological optimizations haven't given satisfactory results because the production cost remains so much expensive. In fact, it seems to be important to take into account physical aspects in order to better understand the mechanical phenomena during the crushing process. Some authors have developed predictive models for physical parameters (friction and wear coefficients). However, these models are based on chemical approaches and are not appropriated to reproduce mechanical effects; force chains distribution in the granular packing and fragmentation modes of grains. Our goal, through out this study, consists to develop a discrete model for granular material behavior subjected to compressive stresses by using DEM [START_REF] Aström | Fragmentation of grains in a two-dimensional packing[END_REF]. In order to better understand the mechanical phenomena, we must understand how the contact force chains are distributed into the granular media subjected to compressive stresses according to the granulometry of the material. In addition, we must study the effect of the granulometry and its impact on the fragmentation modes.

DISCRETE APPOACHE OF MECHANICAL BEHAVIOR

The conventional DEM allows to model really deformable particles as well as complex shapes (from the ellipsoid to the polygon). Here, we have studied the simple case of non-deformable and non penetrable particles in 2D with the computational software MULTICOR, developed by Fortin et al. [START_REF] Fortin | Numerical simulation of granular materials by an improved discrete element method[END_REF]. The coordinates and the rotations of Euler are the configuration parameters q. The gyroscopic and centrifugal terms are equal to zero in 2D. The matrix of generalized mass M of the system doesn't depend on q that is diagonal block. The mechanical equation ( 1) can be written in the following form:

( ) α R t , q , q F q M ext + = & & & (1)
where F ext represents the known external forces and α R the unknown interior forces related to contact reactions and α the contact number. In a system composed of p heterogeneous particles, the critical parameter for the modeling time is the maximum number of interactions between particles. The more the interaction range is important the more we have to test the possible interactions between particles. MULTICOR software uses the partitioning method coupled to a connectivity table [START_REF] Fortin | Selecting contact particles in dynamics granular mechanics systems[END_REF]. The introduction of Coulomb's friction leads to a non linear problem which can't be solved by a linear programming method. Unlike the usual approach, the bipotential method leads to a single variational principle and an inequality [START_REF] Fortin | Numerical simulation of granular materials by an improved discrete element method[END_REF].

STRESS TENSOR FOR GRANULAR MEDIA

The fragmentation mechanisms of grains subjected to compressive loading depend on stresses acting into the granular media. The knowledge of stress field gives the possibility to predict the fragmentation with respect the uniaxial compression strength of the granular media. Modelizing grains as rigid solids in punctual contact for each grains p occupying the volume V p , the expression of mean stress tensor is given by the following equation [START_REF] Fortin | Construction of an average stress tensor for a granular medium[END_REF][START_REF] Guessasma | Numerical modelling of mechanical tests using discrete element method[END_REF].

( )

        ⊗ + - ⊗ = ∫ ∑ p V c p p r q dV q g q V 1 α α α ρ σ & & ( 2
)
where g denote the gravity, q the position vector, q& & its acceleration and α r the contact reaction acting on the particle at the point of contact α c q .

XXII ICTAM, 25-29 August 2008, Adelaide, Australia

The distribution of contact force chains into the granular material resulting from the compressive loading induces stress fields which are responsible of the fracture of grains (figure 1a). To predict the fragmentation of grain, we have considered the Coulomb-Mohr criterion. This criterion states that failure occurs when the two principal stresses are outside the domain of admissible stresses. Several experimental tests using granular materials, under compressive loading, have been carried in order to understand the different fragmentation modes. In our case, we have considered the fragmentation mode of grains which is represented in figure 1b [6].

NUMERICAL SIMULATION OF GRANULAR MEDIA UNDER PRESSURE

We investigate a numerical two-dimensional model for the fragmentation of a granular material under pressure (Figure 2). The calculations are carried out with MULTICOR software. We consider that the behavior of granular material is perfectly rigid with density equal to 2000kg/m 3 and friction coefficient about 0.1. The granular packing contains 32 grains having diameter equal to 15mm and uniaxial compression strength about 5MPa. The granular system is subjected to a pressure velocity about 0.2m/s. The results presented in figure 2 show the fragmentation of granular material and the corresponding contact force chains distribution. The grains subjected to compressive stresses close to the uniaxial compression strength are fragmented following their principal directions. The uniaxial compression stress at the middle of the granular packing plotted with respect the time shows the evolution of stress field according the fragmentation phases.

CONCLUSIONS

We have shown that the contact force chains are crucial to the fragmentation of granular material under pressure. By considering the gravity forces after fragmentation, the rearrangement and distribution of grains seems more realistic. In addition, the overall behavior of granular material reproduces phenomena which are observed in granular media like vault effect or dilatancy.

Figure. 1 .

 1 Figure.1. (a) Simulation of crushing process -(b) Numerical model of fragmentation

Figure. 2 .

 2 Figure.2. (a) Forces chains distributions -(b) Uniaxial compression stress.


Martino Alessandrini

Brecht Heyde

Szymon Cygan

Maxime Sermesant

Hervé Delingette

Olivier Bernard

De Craene

Jan D'hooge

B R E C H Th E Y D E * ,S Z Y M O Nc Y G A N †

Mathieu De Craene

Initial results of a direct comparison in silico based on a new evaluation pipeline

Elastic registration vs. block matching for quantification of cardiac function with 3D ultrasound:

I. INTRODUCTION

Cardiac strain quantification from 3D ultrasound is an active area of research. Commercially available solutions for strain estimation are mainly based on block-matching (BM). Recently, techniques based on elastic registration [START_REF] Heyde | Influence of the grid topology of free-form deformation models on the performance of 3d strain estimation in echocardiography[END_REF] havealso been proven feasible although still confined to the research arena and less extensively tested. Hereto, initial comparison studies recently showed that elastic registration is potentially more accurate BM as far as 2D strain is concerned [START_REF] Heyde | Elastic image registration versus speckle tracking for 2-d myocardial motion estimation: A direct comparison in vivo[END_REF], [START_REF] Jasaityte | Comparison of a new methodology for the assessment of 3d myocardial strain from volumetric ultrasound with 2d speckle tracking[END_REF]. Surprisingly, a thorough and reproducible comparison between the two solutions on 3D ultrasound is to date still missing likely, in part, due to the lack of a solid benchmarking framework.

In this context, we have recently proposed a pipeline for the quality assurance of strain estimation techniques in silico.T he pipeline generates synthetic 3D echocardiographic sequences which looks extremely similar to real ultrasound recordings, while the synthetic cardiac motion is fully controlled by an advanced electro-mechanical model of the heart [START_REF] De Craene | Ultra-realistic synthetic echocardiographic sequences for quality assurance of strain estimation algorithms[END_REF].

In this study, the new pipeline is employed to perform an initial comparison between BM and elastic registration techniques. Hereto, five synthetic datasets were generated including one healthy motion pattern and four ischemic. The two techniques were then compared in terms of tracking accuracy and strain estimation. The capability of discriminating between healthy and akinetic segments were also assessed. The paper proceeds as follows. Sect. II briefly introduces the evaluation pipeline, presents the two algorithms considered and describes the adopted performance metrics. Sect. III reports the results while Sect. IV discusses the limitations of the study, draws the conclusions and presents future perspectives.

II. MATERIALS A N D METHODS

A. Evaluation Pipeline

The evaluation pipeline was recently presented in [START_REF] De Craene | Ultra-realistic synthetic echocardiographic sequences for quality assurance of strain estimation algorithms[END_REF]. It combines state-of-the art solutions in the fields of electromechanical modeling [START_REF] Marchesseau | Preliminary specificity study of the Bestel-Clement-Sorine electromechanical model of the heart using parameter calibration from medical images[END_REF] and ultrasound simulation [START_REF] Gao | A fast convolutionbased methodology to simulate 2-d/3-d cardiac ultrasound images[END_REF] within an original framework that exploits a real ultrasound recording to learn and simulate realistic speckle textures. The simulated images show typical artifacts that make motion tracking in ultrasound challenging. Moreover, the ground-truth displacement field is available voxelwise and is fully controlled by the electromechanical model.

By varying the value of the parameters of the E/M simulator [START_REF] Marchesseau | Preliminary specificity study of the Bestel-Clement-Sorine electromechanical model of the heart using parameter calibration from medical images[END_REF] a library of five cases was created including one healthy heart and four ischemic. The healthy case was generated by assigning normal contractility and stiffness values to all heart segments, as defined by the American Heart Association (AHA) 17 segments model [START_REF] Cerqueira | S tan d ar d ized m y o car d ial s eg m en tatio n an d n o m en clature for tomographic imaging of the heart: A statement for healthcare professionals from the cardiac imaging committee of the council on clinicalc a r d i o l o g y of the american heart association[END_REF]. The ischemic cases were then simulated by modifying peak contractility and stiffness values in diseased segments by differentiating between fully ischemic and mildly ischemic segments. The four ischemic cases corresponded to occlusions of the distal segment of the Left Anterior Descending (LADdist) artery, of the proximal occlusion of the same artery (LADprox), of the Right Coronary Artery (RCA) and of the Left Circumflex artery (LCX), cf. [START_REF] De Craene | 3d strain assessment in ultrasound (straus): A synthetic comparison of five tracking methodologies[END_REF] for more details.

As y n t h e t i cs y s t e mw i t ht h ef o l l o w i n gc h a r a c t e r i s t i c sw a s implemented in the ultrasound simulator [START_REF] Gao | A fast convolutionbased methodology to simulate 2-d/3-d cardiac ultrasound images[END_REF]: sampling frequency 50 MHz, center frequency of the probe 3.3 MHz with a-6d Br e l a t i v eb a n d w i d t ho f6 5 % ,s t a t i cf o c u s i n ga t8 0m m when transmitting and dynamically focusing on receive. The simulated images consisted of 107×80 lines in azimuth and elevation direction over an angle of 76×76 degrees, resulting in a frame rate of 30 Hz assuming parallel beam forming and ECG gating. After scan-conversion each volume measured 208×176×208 pixels 3 (voxel size 0.6×0.8×0.5 mm 3 ). An example of simulated sequence is provided in Fig. 1. An animation showing one full sequence can be found at http://bit.ly/sim us.

B. Considered Algorithms

All five simulated sequences were processed by the two motion estimation techniques considered.

1) Elastic Registration:

The technique in [START_REF] Heyde | Influence of the grid topology of free-form deformation models on the performance of 3d strain estimation in echocardiography[END_REF] was considered as representative of the elastic registration family. In its formalism, the inter-frame displacement field u = u(x, y, z) is expressed as:

u = i,j,k µ i,j,k β x -κ i x σ x β y -κ j y σ y β z -κ k z σ z (1)
being β the third order B-spline, κ ξ the node spacing and σ ξ the node span. Parameters µ i,j,k are optimized with the LBFGSB optimizer by minimizing the cost function:

E = r ∥I n -I n+1 (T(r))∥ 2 + α r ∂ 2 T(r) ∂r 2
2

(2)

where T(r)=r + u(r) and α is the hyperparameter to modulate the influence of the bending energy.

As opposed to previous works where control points were put on a Cartesian grid, an original anatomical grid topology is employed in [START_REF] Heyde | Influence of the grid topology of free-form deformation models on the performance of 3d strain estimation in echocardiography[END_REF], specifically adapted to the geometry of the considered ventricular shape. With respect to the Cartesian topology, the novel technique was shown to be competitive in terms of estimation accuracy while computationally more effective, cf. [START_REF] Heyde | Influence of the grid topology of free-form deformation models on the performance of 3d strain estimation in echocardiography[END_REF].

2) Block matching: Block matching employed 3D Normalized Cross-Correlation (NCX) as similarity metric. Prior to localizing the maximum of the resulting NCX it was interpolated in 3D by recursively (3 times) interleaving interpolates between every element using the cubic B-spline method, which resulted in 8-times improved spatial resolution of displacements estimates. The search range for each block matching was set globally to cover maximum physiological velocities in cardiac motion of 12cm/s. Dimensions of matched blocks (kernels) were chosen experimentally by assessing the average and median errors of displacements estimation for various kernel sizes. Here the kernel was 14.3 × 18 × 12.3 mm 3 .N os moot hi ng,re gul a ri z a t i onorqua l i t yc he c kme a s ure s were implemented.

C. Performance Evaluation

For each sequence, the result of each estimator was applied to progressively track the 5387 nodes of the first simulation mesh over the whole cardiac cycle. Tracking accuracy was then evaluated by measuring the Euclidean distance to the nodes of the ground truth meshes. In particular, for each technique we evaluated global tracking error (i.e. inclusive of all nodes and time instants), instantaneous tracking error (i.e. inclusive of all nodes at a specific time instant) and local tracking error maps evaluated at end-systole (ES).

From the sets of volumetric meshes, strain values were obtained by approximating the displacement within each tetrahedron as an affine transform F.T h eG r e e n -L a g r a n g es t r a i n tensor E =1/2 F T F -I was then computed and projected on the local cardiac coordinate system of longitudinal (Lstrain), radial (R-strain) and circumferential (C-strain)d i r For each algorithm we evaluated the accuracy in recovering ES strain (the most relevant from a clinical perspective). Hereto, Pearson correlation coefficient together with the bias and standard deviation returned by the Bland-Altman (BA) analysis were used. For each correlation value the p-value was computed testing the hypothesis of no correlation. The statistical significance of each reported bias was measured with a t-test. Fishers z-transform was used to compare the strengths of different correlations. A t-test was used to compare the biases returned by the BA analysis. Differences were considered statistically significant if p<0.05.

III. RESULTS

A. Tracking accuracy

Boxplots of the global tracking error are reported in Fig. 2. Elastic registration was on average more accurate than block matching on all sequences considered. On the entire dataset, global tracking error was (µ ± σ)1 . 1 5 3 ±0.625 mm for elastic registration and 1.612±1.161 for block matching. The same behavior was observed when considering the time evolution of the tracking error, cf. Fig. 3.

Besides smaller average errors, elastic registration returned more robust estimates (i.e. with a higher snr = µ/σ). This is clearly due to the intrinsic spatial regularization of the estimated displacement fields. As a consequence, smoother error maps were obtained with elastic registration, cf. Fig. 4 .

B. Strain accuracy

Correlation plots and Bland-Altman plots of the computed ES circumferential strain are reported in Fig. 5. All sequences were considered at once, i.e. 85 sample points (17 segments times 5 sequences) were included in the analysis for each algorithm. The values found for all the strain components are summarized in Table I. For all the strain components elastic registration returned higher correlation coefficients, lower bias and smaller limits of agreement as compared to block matching. Nevertheless, observed differences were not always statistically significant (cf. Table I). This is likely due to the limited number of samples included in this initial study. For both algorithms, the most reliable strain component was the circumferential one, while significant biases and low correlation values were obtained for the radial component. In agreement with previous studies, this confirms that obtaining reliable estimates of radial displacements is still an open issue.

As a measure of the sensitivity of the two techniques to the presence of akinetic regions, we evaluated how effectively computed strain values could cluster healthy and diseased segments. The results are summarized in Fig. 6. We used the Mahalanobis M =(µ 1 -µ 2 )/(σ 1 +σ 2 ) to evaluate the distance between the two clusters. For all strain components elastic registration produced the best inter-cluster separation (cf. Table II). Nonetheless, both methods performed poorly when Rstrain was considered. In line with previous studies, this result confirms that accurate evaluation of radial displacements with ultrasound is still an open issue. The paper presented an initial comparison between elastic registration and block matching techniques for cardiac deformation and strain estimation with 3D ultrasound making use of a novel evaluation pipeline recently proposed by the (f) R-strain Fig. 6. In the x-axis the reference strain for healthy (in black) and ischemic (in red) segments. In the y-axis the values returned by the estimator. The two best-fit Gaussian curves are reported as an indicator of the separation between the two groups.

authors. The main limitation of this study was that an in-house implementation of block matching was employed where i) no spatial nor temporal smoothing were included, ii) computational efficiency was highly sub-optimal (≈ 6hp e rs e q u e n c e ) . Hereto, in the future it would be of interest to involve commercially available implementations of speckle tracking int h e comparison. Moreover, given the reduced size of the dataset, obtained differences were not always statistically significant. Am o r ee x t e n s i v ee v a l u a t i o ni st h e r e f o r en e e d e d ,i n c l u d i n ga wider range of motion patterns and image qualities.

Despite these limitations, this study brings further evidence that elastic registration might represent a valuable alternative to speckle tracking for a reliable estimation of cardiac dynamics and that this novel technique is hence ready for more extensive clinical trials. Future studies include employing the proposed benchmarking pipeline to contrast a wider range of state-of-the art motion estimation solutions.

Fig. 1 .

 1 Fig. 1. End diastolic and end-systolic frames from the healthy simulation. In red the boundary of the ground-truth meshes providing reference values of displacement and strain.

  e ctions according to d = d T Ed.B e i n gd =[ d 1 ,d 2 ,d 3 ] T the unit vector for the anatomical coordinate of interest. Segmental strain values were then obtained by averaging the point estimates for each AHA segment.

Fig. 2 .

 2 Fig. 2. Global tracking errors for elastic registration (in red) and block matching (in blue) on the five considered sequences. Markers in each box correspond to 25-th, 50-th and 75-th percentile values. Values outside the range +/-3σ are classified as outliers (red pluses).

Fig. 3 .

 3 Fig. 3. Tracking error over time. Bold time curves report median errors. Ver tical s eg m en t d en o te 2 5 -th an d 7 5 -th p er cen tile valu es .

Fig. 4 .

 4 Fig. 4. Error maps for elastic registration (a) and block matching (b) at endsystole of the healthy sequence. Error values are displayed on the ground-truth geometry.

Fig. 5 .

 5 Fig. 5. Correlation (a,b) and Bland-Altman analysis (c,d) for circumferential ES strain estimated with the two techniques considered. The correlation plot also reports the expression of the line fitted in the least-squares sense.

TABLE II

 II

	.	M AHALANOBIS DISTANCE BETWEEN CLUSTERS
	CORRESPONDING TO HEALTHY AND DISEASED SEGMENTS.
			Reference	BM	Registration
		C-strain	1.64	1.09	1.29
		R-strain	1.54	0.46	0.74
		L-strain	1.12	0.69	0.87
	IV. DISCUSSION AND CONCLUSION

TABLE I .

 I C ORRELATION COEFFICIENT,BLAND-ALTMAN BI AS µ AND BLAND-ALTMAN LI MI TS OF AGREEMENT σ IN THE ESTIMATION OF ES STRAIN. THE p-VA L U E B E T W E E N B R AC K E T S R E P O RT S T H E S TAT I S T I C A L S I G N I FI C A N CE OF THE REPORTED VALUE.THE SYMBOL† DENOTES STATISTICAL DIFFERENCE (p<0.05) BETWEEN THE TWO ALGORITHMS.

			R-strain			C-strain			L-strain
		ρ (p-value)	µ (p-value)	σ	ρ (p-value)	µ (p-value)	σ	ρ (p-value)	µ (p-value)	σ
	ELASTIC REGISTRATION	0.74 (<0.001)	-2.55 (0.01) †	17.58	0.94 (<0.001) †	0.98 (<0.001) †	4.94	0.71 (<0.001)	-0.28 (0.25) †	4.44
	BLOCK-MATCHING	0.56 (<0.001)	13.12 (<0.001) †	21.67	0.77 (<0.001) †	-1.12 (0.03) †	9.30	0.61 (<0.001)	-0.95 (0.002) †	5.29

ACKNOWLEDGMENT

The work of Martino Alessandrini was funded by the Research Foundation Flanders (FWO).