domain: spi phys

text:
Mohamed-Amine Chadil
email: achadil@imft.fr

Stéphane Vincent
email: stephane.vincent@u-pem.fr

Jean-Luc Estivalèzes
email: jean-luc.estivalezes@onera.fr

M.-A Chadil

J.-L Estivalèzes

Accurate estimate of drag forces using particle-resolved direct numerical simulations

Keywords: Accurate drag force calculation, Finite size particles, Immersed boundary techniques, Penalty methods, Fixed arrangements of spheres

An accurate force estimate for finite-size particle simulations is proposed based on Lagrange extrapolation of third order, coupled with a Taylor interpolation of same order, to estimate pressure and viscous constraints on the surface of particles. The main point of our approach is to upwind the interpolation support in the normal direction to the fluid/solid interface so as to use only fluid values to estimate forces. Also, detailed validations of forces are considered for estimating accuracy and convergence order of the method on various incompressible motions such as the flow around an isolated particle at various Reynolds numbers and flows across packed spheres under Faced-Centered Cubic, random and bi-disperse arrangements.

Introduction

In the framework of finite-size particle motions, the numerical simulation of a particulate flow interacting with a surrounding fluid can be investigated following two different numerical strategies: unstructured or structured grids. This important choice is in particular motivated by the instantaneous description of the evolving complex shape represented by the interface between the carrier fluid and a set of moving particles. On the one hand, the more natural discretization seems to be the implementation of an unstructured body-fitted grid to simulate the fluid area in the two-phase particle flow [START_REF] Hu | Direct Simulation of Fluid Particle Motions[END_REF] [2] [3] [START_REF] Cho | Direct numerical simulation of fluid flow laden with many particles[END_REF]. Building such a finite-volume or finite-element mesh in three-dimensions is not easy and requires automatic remeshing as the solid particles move according to time and space. The remeshing process at each calculation step is time consuming and can be very difficult to manage automatically in computer softwares when the global shape of the fluid-solid interface is changing at each calculation step [START_REF] Dobrzynski | Anisotropic Delaunay Mesh Adaptation for Unsteady Simulations[END_REF]. On the other hand, it can be imagined to develop a fixed structured grid approach to simulate particulate flows. With this method, the mesh is not adapted to the fluid-solid interfaces and includes both phases. On a mesh point of view, this approach is simpler than the previous one. The difficulty lies in taking into account the presence of particles in the fluid whose interface is not explicitly tracked by the mesh that does not conform to the fluid/particle interface. This type of modeling and numerical problem belongs to the class of fictitious domains [START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow[END_REF] [START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF]. The penalty modeling strategy [START_REF] Vincent | An adaptative augmented Lagrangian method for three-dimensional multi-material flows[END_REF] [9] [START_REF] Vincent | A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows[END_REF] developed hereafter is based on this approach and will be reported in the next section.

One major interest of finite-size particle motion simulation is to provide a local estimate of all flow characteristics (velocity, pressure, viscous stress) together with a local description of the particle-fluid interface. The resolved scale particle motion does not require any force modeling nor interaction model as soon as the mesh is refined enough. In the present work, we will demonstrate that having at least 5 points in the boundary layer attached to a given particle allows to recover all the physics of the particle flow without using any force or interaction model (drag, lift, lubrication, etc). As a consequence, the finite-size particle approach can be considered as a kind of fully resolved Direct Numerical Simulation (DNS) of the particulate flow as soon as no modeling is required to solve the problem. The mesoscopic or macrocopic particle flow models (Eulerian-Lagrangian, Eulerian-Eulerian) do require the knowledge of interaction forces to close specific particle-motion interaction terms. For an isolated particle in an infinite medium, the drag force law of Schiller and Naumann [START_REF] Schiller | Über die grundlegenden Berechnungen bei der Schwerkraftaufbereitung[END_REF] is well known and often used in large scale models. For fixed and moving beds of particles, we can cite the correlations based on experiments proposed by Ergun [START_REF] Sabri Ergun | Fluid flow through packed columns[END_REF] and Wen&Yu [START_REF] Wen | Mechanics of fluidization[END_REF] and also Gobin [14] who proposed a correlation based on these two correlations. As soon as the solid fraction is high and the particle size or shape of particles is not constant (bi-dispersed flows, spheroidal particles, etc), drag and lift laws have to be designed. This is the main objective of the present work, i.e. providing a DNS framework with accurate force calculation in order to finally obtain new force laws for various particulate motions.

Numerous numerical works have been devoted to performing DNS of finitesize particulate flows and obtaining resulting drag force laws from the macroscopic analysis of local flow motions in the vicinity of the particle surface. The first class of numerical models, generally investigated for fixed arrangements of particles, is the Lattice Boltzmann approach. It was used in numerous works by Ladd [START_REF] Anthony | Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation[END_REF] [START_REF] Anthony | Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 2. Numerical results[END_REF] for particle suspensions, by Hill et al. [START_REF] Reghan | Moderate-Reynolds-number flows in ordered and random arrays of spheres[END_REF] for monodiperse face cubic centered (FCC) array of spheres, by Hoef et al. [START_REF] Van Der Hoef | Lattice-Boltzmann simulations of low-Reynolds-number flow past mono-and bidisperse arrays of spheres: results for the permeability and drag force[END_REF] and beetstra et al. [START_REF] Beetstra | Drag force of intermediate Reynolds number flow past mono-and bidisperse arrays of spheres[END_REF] for random monodisperse and bidisperse array of spheres to cite a few. Another class of very popular methods for handling finite-size particle is the Immersed Boundary Method (IBM) coupled with incompressible Navier-Stokes equations initially developped for particulate flows by [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF] [START_REF] Uhlmann | Interface-resolved direct numerical simulation of vertical particulate channel flow in the turbulent regime[END_REF]. For random monodisperse array of spheres, the work of Tenneti et al. [START_REF] Tenneti | Drag law for monodisperse gas-solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres[END_REF] is very interesting as it covers a wide range of solid fractions α d = 0.1 -0.4, and also Reynolds number up to 300. At the end, a new correlation is proposed for drag force laws. The last interesting class of numerical approaches is the body fitted mesh method, that is restricted to fixed array of particles as it is impossible to generate automatically a three-dimensional mesh that adapts to the motion of particles. Among the most interesting works in the field, we can cite the simulations and analysis of Massol et al. [START_REF] Massol | Simulations numériques d'écoulements travers des réseaux fixes de spheres monodisperses et bidisperses, pour des nombres de Reynolds modérés[END_REF] for monodispersed FCC arrays of spheres.

In all these numerical approaches, the drag force can be deduced in different ways:

• with Darcy penalty methods [START_REF] Nirina | Etude numérique des interactions hydrodynamiques fluides/solides : application aux lits fluidisés[END_REF], it can be obtained directly from the source term used to accelerate the fluid ouside the particles. • with the LBM, the drag force exerted by the fluid on the particles is calculated according to the momentum exchange algorithm of Ladd [START_REF] Anthony | Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation[END_REF]. • In the IBM approach, the forces exerted by the fluid on the particles can be deduced from the reaction IBM force imposed in the momentum equations to satisfy a solid behaviour [START_REF] Breugem | A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows[END_REF]. • with the Volume Penalty Method (VPM), Bizid [START_REF] Bizid | Penalty methods for turbulent flows interacting with obstacles[END_REF] uses Taylor extrapolations to get the pressure and viscous stress projection on particle surface.

In the framework of macroscopic Eulerian-Eulerian or mesoscopic Eulerian-Lagrangian models, the previously cited drag force laws are used in many industrial processes such as fluid catalytic cracking reactors [START_REF] Ali | Dynamic modeling and simulation of a risertype fluid catalytic cracking unit[END_REF], gas phase polymerisation reactors [START_REF] Gobin | Fluid dynamic numerical simulation of a gas phase polymerization reactor[END_REF], chemical looping [START_REF] Zhang | Modelling and simulation of chemical looping combustion process in a double loop circulating fluidized bed reactor[END_REF] or fluidized beds [START_REF] Vincent | A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows[END_REF] [START_REF] Ozel | Particle resolved direct numerical simulation of a liquidsolid fluidized bed: Comparison with experimental data[END_REF]. However, convergence studies and accuracy analysis of force calculation have rarely been considered.

Our major contribution is to propose an accurate force estimate for finitesize particle simulations. Even if all full DNS of particle flows extract forces on particles [START_REF] Vincent | A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows[END_REF] [20] [START_REF] Randrianarivelo | A DNS approach dedicated to the analysis of fluidized beds[END_REF] [2], few works report on how practically these forces can be calculated and what is the accuracy or convergence order of the forces. Among the wide literature devoted to full DNS of particle motions, we can cite the work of Bizid [26] [31] who uses Taylor extrapolation to estimate pressure and viscous constraints on the particle surface. This approach is of low accuracy as it utilizes velocity and pressure values that can be inside the particle volume. The most advanced work on force calculation on immersed interfaces is due to Zastawny et al. [START_REF] Zastawny | Derivation of drag and lift force and torque coefficients for non-spherical particles in flows[END_REF] with an improved mirroring immersed boundary method. In the present work, a new force calculation is proposed based on Lagrange extrapolation of third order, coupled with a Taylor interpolation of same order. The main point of our approach is to upwind the interpolation support in the normal direction to the fluid/solid interface so as to use only fluid values to estimate forces. Detailed validations of forces are considered for estimating accuracy and convergence order of the proposed method.

The article is structured as follows. A presentation of the models and numerical methods is first proposed in section 2, paying attention to describe fictitious domain and penalty methods used to model and approximate incompressible particulate motions. In the third section, a new force calculation for immersed interfaces is proposed, with discussions conducted on order of approximations and associated accuracy. Validations for flows interacting with isolated spherical particles at various Reynolds number are presented in the fourth section. Section 5 is devoted to simulations and validations of flows through fixed arrangements of mono-and bi-dispersed spheres. Finally, conclusions and perspectives are drawn.

Model and numerical methods

Fictitious domain approach

The modeling and simulation of moving objects (bubbles, droplets, solid particles) interacting with a carrier fluid is impossible to realize with unstructured meshes as soon as these objects deform or move in a 3D geometry. The commonly developed alternative approach consists in simulating this kind of flow on a fixed grid and to locate the interface thanks to an auxiliary phase function such as Volume Of Fluid or Level Set functions [START_REF] Sarthou | A second-order curvilinear to cartesian transformation of immersed interfaces and boundaries. application to fictitious domains and multiphase flows[END_REF]. The concept that disconnects the interface motion and the mesh used to solve the conservation equations is called fictitious domain approach [START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF] [START_REF] Randrianarivelo | Numerical modelling of the solid particle motion using a new penalty method[END_REF]. Indeed, from the motion equation point of view, the interface is not known, only its presence is taken into account thanks to a volume auxiliary function. In these approaches, the interface tracking and the associated building of the phase function is of primary importance.

As proposed in [START_REF] Vincent | A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows[END_REF], incompressible two-phase flows involving a carrier fluid and a solid phase can be modeled on a fixed mesh with fictitious domain approaches by considering the incompressible Navier-Stokes equations together with a phase function C describing the particle phase shape evolutions through an advection equation on the corresponding phase function. As explained by Kataoka [START_REF] Kataoka | Local instant formulation of two-phase flow[END_REF] for fluid/fluid two-phase flows and Vincent [START_REF] Vincent | A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows[END_REF] for particle flows, the resulting model takes implicitly into account the coupling between different phases separated by resolved interfaces, i.e. larger than the mesh cell size. The motion equations reads

∇ • u = 0 (1) ρ ∂u ∂t + (u • ∇) u = -∇p + ρg + ∇ • µ(∇u + ∇ t u) + F si + F m (2) ∂C ∂t + u • ∇C = 0 (3)
where u is the velocity, p the pressure, t the time, g the gravity vector, ρ and µ respectively the density and the viscosity of the equivalent fluid. The twoway coupling between particle and fluid motions is ensured in the momentum equations by the presence of a solid interaction force F si [START_REF] Brändle De Motta | Fictitious Domain Methods and Penalty Techniques for the Simulation of Turbulent Particulate Flows[END_REF] [37] which is not considered in the present work as only fixed particles are dealt with. The source term F m is used to impose a flow rate to the fluid.

The one-fluid model is almost identical to the classical incompressible Navier-Stokes equations, except that

• the local properties of the equivalent fluid (ρ and µ) depend on C. In the present work, an arithmetic average is used for density and an harmonic average is considered for viscosity [START_REF] Vincent | A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows[END_REF]. • the interface localization requires the solving of an additional equation on C. Instead of solving this equation on the Eulerian mesh, which is source of numerical diffusion or tearing of interfaces, a Lagrangian representation is preferred. A specific mesh (linear elements in 2D and set of triangles in 3D) is considered for the particle surface S. Equation ( 3) is reformulated as dX b dt = V p , with X b the centroid of the spherical particle and V p the velocity of the particle interpolated with surrounding Eulerian velocities coming from the solving of equations (1-2). The approximation of the Lagrangian tracking of X b is detailed in [START_REF] Vincent | A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows[END_REF]. Once X b is known, the position of each particle surface mesh element is also known. The phase function C is automatically build by projecting the Lagrangian particle mesh onto the Eulerian mesh [START_REF] Sarthou | A second-order curvilinear to cartesian transformation of immersed interfaces and boundaries. application to fictitious domains and multiphase flows[END_REF]. For non-spherical particles the rotational motion has to be considered. It is for example solved with Quaternions [START_REF] Eberly | Quaternion algebra and calculus[END_REF]. In the present work, only fixed arrangements of particles are dealt with. As a consequence, the Lagrangian particle tracking is not considered even if the Lagrangian mesh for representing particles is used to calculate forces.

Satisfying the incompressible and solid constraints in fluid and particles requires developing a specific model. Two penalty approaches are proposed and detailed in the next section to tackle with these constraints.

Penalty methods

As explained in the previous section, the one-fluid model and the fictitious domain approach formulated for dealing with particle flows require to consider each different phase (fluid, solid) as a fluid domain with specific material properties (density and viscosity for an isothermal flow). Each sub-domain is located by a phase function C.

In addition to local physical characteristics of the fluid that change over time due to particle motion, the local constraints that have to be fulfilled change potentially at each time step. Indeed, a given point or cell can be inside a fluid zone at time n∆t and can be solid at the next time step (n+1)∆t. Here, n is the time index and ∆t the associated time step. Two different numerical methods are used to satisfy in a coupled way and at the same time the fluid incompressibility and the solid behavior:

• Ensuring the solid behavior in the solid zones where C = 1 requires to define a specific rheological law for the rigid fluid part without imposing the velocity, as the particle velocities are not always known a priori in particulate motions (particle sedimentation, fluidized beds, turbulence particle interaction). A specific model is implemented for handling the solid particle behavior in the one-fluid Navier-Stokes equations. It is based on a decomposition of the viscous stress tensor and on a penalty method that acts on the viscosity which tends to large values in the particles [START_REF] Caltagirone | Tensorial penalisation method for solving Navier-Stokes equations[END_REF] to implicitly impose the solid behavior and also the coupling between fluid and solid motion. For fixed particles, the velocity of the cell containing the centroid of the particle is imposed equal to zero. The viscous penalty method propagates the zero velocity in the whole solid medium. • Following a similar walkthrough as in the work on Stokes and Navier-Stokes equations proposed by Fortin and Glowinski [START_REF] Fortin | Méthodes de lagrangien augmenté. Application à la résolution numérique de problèmes aux limites[END_REF], an augmented Lagrangian method is applied to the unsteady Navier-Stokes equations dedicated to particulate flows. The main objective is to deal with the coupling between velocity and pressure and to satisfy the fluid and solid constraints at the same time. Starting with u * ,0 = u n and p * ,0 = p n , the augmented Lagrangian solution reads while

||∇ • u * ,m || > , solve (u * ,0 , p * ,0 ) = (u n , p n ) ρ u * ,m -u * ,0 ∆t + u * ,m-1 • ∇u * ,m -∇(r∇ • u * ,m ) = -∇p * ,m-1 + ρg + ∇ • [µ(∇u * ,m + ∇ T u * ,m )] + F si + F m p * ,m = p * ,m-1 -r∇ • u * ,m (4)
where r is an augmented Lagrangian parameter used to impose the incompressibility constraint, m is an iterative convergence index and a numerical threshold controlling the constraint. The augmented Lagrangian method is a kind of penalty technique: if r → +∞, the incompressibility is imposed but the solving of the linear system is difficult with iterative solvers while r → 0 does not act on the fluid constraint and keeps the conditioning of the matrix unchanged. Usually, a constant value of r is used, for example equal to the average between the minimum and maximum eigenvalues of the linear system for Stokes flows [START_REF] Fortin | Méthodes de lagrangien augmenté. Application à la résolution numérique de problèmes aux limites[END_REF]. From numerical experiments, optimal values are found to be of the order of ρ i and µ i in each phase (fluid or solid) to accurately solve the motion equations in the related zone [8] [9]. Algebraic improvements have also been proposed by Vincent [START_REF] Vincent | Augmented Lagrangian and penalty methods for the simulation of two-phase flows interacting with moving solids. Application to hydroplaning flows interacting with real tire tread patterns[END_REF] to automatically estimate the local values of r. In the present work, a constant value of r will be used. Its magnitude will be discussed in the validation sections.

Discretization schemes and solvers

The schemes and solvers used in the present work are detailed in [START_REF] Vincent | A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows[END_REF]. The mass and momentum conservation equations, containing the viscous and augmented Lagrangian penalty terms, are discretized with implicit finite volumes on structured staggered meshes. The time derivative is approximated with a first or second order Euler scheme while the inertial, viscous and augmented Lagrangian terms are discretized with a second-order centered scheme. All fluxes are written at time (n + 1)∆t, except the non-linear inertial term that is linearized at first or second order as follows

• u • ∇u ≈ u n • ∇u n+1 for first order linearization • u • ∇u ≈ 2u n -u n-1 • ∇u n+1
for second order Adams-Bashforth like linearization

The obtained linear system is solved with a BiCGSTAB II iterative solver [START_REF] Gustafsson | On first-and second-order symmetric factorisation methods for the solution of elliptic difference equations[END_REF], preconditionned with a Modified and Incomplete LU approach [START_REF] Van Der | Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of non-symmetric systems[END_REF]. All the code is working on massively parallel computers by using MPI devices and exchanges [START_REF] Vincent | A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows[END_REF].

Lagrangian extrapolation of forces for immersed boundary methods

The drag force due to the fluid-solid interaction at a surface S of the solid phase is:

F D = F p + F v (5)
where the pressure force F p and the viscous force F v are:

F p = S -p n dS (6)
F v = S 2µσ.n dS (7)
𝑛 𝑙 The computation of these forces consists in discretizing S on a set of N elements (triangles in 3d and segments in 2d) called Lagrangian mesh (see figure 1), such that:

Fluid Solid

F p ≈ N l -p l n l dS l (8)
F v ≈ N l 2µσ l .n l dS l (9)
where p l , σ l and n l are respectively the pressure, strain tensor and outgoing normal vector at the center C l of the l th element of the Lagrangian mesh as illustrated in figure 1. Normal n l and element surface dS l are deduced from the coordinates of the nodes constituting the l th element.

Low order naive approach

Given that the pressure field p and the viscous stress tensor σ, used in the drag force computation, are known on the Eulerian mesh and not on the Lagrangian mesh as explained above, the naive approach to overcome this problem consist in interpolating them from the Eulerian mesh on the Lagrangian mesh using a second order Taylor interpolation detailed in Appendix 1.

To validate this approach, we compute the drag force exerted by a uniform Stokes flow (Re = 10 -3 ) on an isolated cylinder (2D) and an isolated sphere (3D), as detailed below.

Uniform Stokes flow past a cylinder

According to [START_REF] Bizid | Development of penalty methods for the simulation of turbulent flows around obstacles[END_REF], a uniform Stokes flow past a cylinder of diameter d = 2m, with the undisturbed velocity being noted U ∞ = 1m/s, is solution of the Brinkman equation -∇p

+ µ∆u i - µ K u i = 0.
In the polar coordinate frame (r, θ) centered on the particle, it reads:

u (r, θ) =    1 r -1 + 2K1(λ) λK0(λ) 1 r + r + 2 λK0(λ) K 1 (λr) cos θ -1 + 1 + 2K1(λ) λK0(λ) 1 r 2 -2 K0(λ) K 0 (λr) + K1(λr) λr sin θ (10) p (r, θ) = µλ 2 -1 + 2K 1 (λ) λK 0 (λ) 1 r -r cos θ (11)
The corresponding drag force is:

F D = 2π Re 2λ 2 + 4 λK 1 (λ) K 0 (λ) e x (12)
where λ =

d 2 4K
is the dimensionless permeability of the porous medium in Brinkman sens, K is the permeability of the porous cylinder, K 0 and K 1 are the modified Bessel functions of rank 0 and 1. For K → 0, the porous cylinder can be likened to an impermeable solid particle.

Uniform Stokes flow past a sphere

For the uniform flow past a sphere in Stokes regime, the velocity and pressure are [START_REF] Keith | An introduction to fluid dynamics[END_REF] :

u (x, y, z) =          U ∞ -3 4 ( d 2 )U ∞ x 2 r 3 + 1 r + 3 4 ( d 2 ) 3 U ∞ x 2 r 5 -1 3r 3 3 4 U ∞ ( d 2 ) 3 r 5 - ( d 2 ) r 3 xy 3 4 U ∞ ( d 2 ) 3 r 5 - ( d 2 ) r 3 xz (13) p = - 3 4 µdU ∞ x r 3 (14)
where

r = x 2 + y 2 + z 2 .
The drag force is analytically given by: Fig. 2: Drag force relative error (%) for the Stokes flow past a cylinder (left), and the Stokes flow past a sphere (right) in terms of Lagrangian mesh refinement -the force is computed with the naive 2 nd order Taylor interpolation method. The size of the Lagrangian particle surface mesh element is l while the size of the Cartesian Navier-Stokes mesh is ∆x.

F D = 3πµdU ∞ e x ( 15
)

Simulations setup

The computational domain used to simulate a uniform Stokes flow past a cylinder (resp. a sphere ) is a square (resp. a cube) of a Length L = 2d, and the spatial discretization, using a regular Cartesian grid called Eulerian mesh, is represented by the number of grid cells across the diameter of the particle

D m = d ∆x = 20.
The exact solution [START_REF] Vincent | A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows[END_REF], [START_REF] Schiller | Über die grundlegenden Berechnungen bei der Schwerkraftaufbereitung[END_REF] (resp. ( 13), ( 14)) for a Stokes flow past a cylinder (resp. a sphere) was implemented at boundary conditions to be able to simulate such a flow in a numerically small domain not extending to infinity as Stokes flow would require. Practically, the considered simulation domain is two particle diameter long in each space direction.

Results

Figure 2 shows that the computation of the drag force with a naive approach is not accurate, the error being about at least 40% even if the surface Lagrangian mesh is refined. This is due to the use of pressure and strain tensor values in the cells containing the fluid-solid interface where the error with the analytic solution is the highest, as illustrated in figure 3 for the uniform flow past a cylinder (2D) and a sphere (3D). This error is due to the fact that the physical characteristics (µ, ρ) in the mixed fluid-solid cells are the average between those of the fluid and those of the solid.

To prevent the use of wrong pressure and velocities in the cells cut by the fluid-particle interface, we have to extrapolate the pressure and the strain tensor from the fluid area far from the particle to the interface. To do so, a Lagrange extrapolation is considered, as detailed below.

f (C l ) = k i=1 f (P i )L i (C l ) , where L i (C l ) = k j =i |C l -P j | |P i -P j | (16)
The k th order of the drag force computation consists now of interpolating (using the k th order Taylor interpolation detailed in Appendix 1) the pressure and the strain tensor on k P i fictitious points built along the normal (see figure 4 for third order), and then using the k th order Lagrange extrapolation [START_REF] Anthony | Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 2. Numerical results[END_REF] to compute them at C l . With ( 5), ( 8), ( 9), the friction force exerted by the fluid on the particle can be calculated.

In figure 4, points P 1 , P 2 and P 3 are the fictitious points used in third order Lagrange extrapolation of the pressure and viscous components of the force on the particle surface. The zone between the particle surface S (solid line) and the dash line corresponds to the Eulerian discrete points that have an influence on the force calculation. Parameter δ is the minimum distance between the first fluid point and S that is required to obtain a correct force estimate.

Note that the order of the drag force computation may not be the same on all the Lagrangian mesh element. Indeed the choice of the computation order located at the centre C l of each Lagrangian mesh element is governed by the fact that all the points used in the Lagrange extrapolation and the Taylor interpolation have to be fully fluid, i.e the fluid-phase indicator function C = 0 on these points.

Considering that the pressure and strain tensor have to be extrapolated from the nearest fluid region in the vicinity of the particle, where the physical characteristics belong to the fluid, to the interface region where the drag force is computed, the first step consists in increasing the distance δ (see figure 4) until the computation of the force is accurate enough. As illustrated in figure 5, the minimum error is reached for δ = ∆x and remains stable for the third and fourth order Lagrange interpolation, for both simulations, namely the Stokes flow past a cylinder (left), and past a sphere (right) . As a consequence, the third order force computation with δ = ∆x will be considered thereafter as it requires less computational effort. The errors obtained with the new high order force calculation method are reported in figure 6. Compared to the naive approach, the error levels for both cylinder and sphere cases are now always less or around 2% as soon as the local size of the Lagrangian surface mesh elements is comparable to ∆x. Thanks to the novel force estimate, the error has been reduced by a factor of 20. In the rest of the work, the value of δ will be chosen to be ∆x if not specified. For validation purpose, the force acting on a particle and in particular the corresponding drag coefficient is recalled. It is defined as [START_REF] Clift | Bubbles, Drops, and Particles[END_REF]:

C d = |F D | 1 2 ρ Ũ 2 A p (17)
where Ũ = |U ∞ e x -U p | is the relative velocity between the particle and the fluid velocity at infinity, U ∞ e x is fluid velocity in the mean flow direction far from the particle, U p is the particle velocity, ν is the fluid kinematic viscosity and A p = π 4 d 2 the cross-sectional area of the particle.

The drag coefficient is dependent on the flow regime determined by the Reynolds number:

Re = Ũ d ν (18)
The drag coefficient in Stokes regime, i.e. when Re → 0, is:

C d = 24 Re (19)
The correlation of drag coefficient for a finite Reynolds number, proposed by Schiller & Naumann [START_REF] Schiller | Über die grundlegenden Berechnungen bei der Schwerkraftaufbereitung[END_REF] is:

C d = 24 Re 1 + 0.15 Re 0.687 (20)
𝑈 ∞ Fig. 7: Illustration of instantaneous stream lines obtained for a uniform flow past an isolated sphere at Re=290. The upstream unperturbed velocity U ∞ is imposed at the left boundary condition.

Simulations setup

The computational domain is first chosen, its lengths being L x = 16d and L y = L z = 8d in each Cartesian direction. The Eulerian mesh refinement is constant in a box of extension [(2d, 3d, 3d) ; (6d, 5d, 5d)] centered around the particle position. Outside this box, the Eulerian mesh is exponentially coarsen from the box to the boundaries of the simulation domain. Accurate drag force calculation needs properly resolved boundary layers around the particle. It is expected that a 5 cell resolution in the boundary layer thickness will be enough at least for the range of Reynolds numbers studied in the present article. This choice is also consistent with the compact support of cells needed to get third order Lagrange extrapolation accuracy (see figure 24 in appendix 1). This gives d √ Re = 5∆x according to scaling laws for laminar dynamic boundary layers.

Therefore, in the box surrounding the particle, the minimum cell size is:

∆x = d 5 √ Re (21)
The inlet boundary conditions is u = U ∞ e x (see figure 7) and Neumann conditions are applied elsewhere.

Study of numerical parameters for the Lagrange extrapolation

As an extension of the work carried out for the Stokes regime, we have performed the same kind of study for a higher finite Reynolds number Re = 100 when inertial effects are important. The effect on the drag force computation of the distance δ between the particle interface points C l and the first Eulerian mesh extrapolation point P 1 has been considered. The role played by the Lagrangian mesh refinement have also been analyzed.

One can observe in figure 8 a similar behaviour of error as already found for Stokes flows (the reader can refer to figure 5 right) indicating that in this range, Reynolds number have little influence on drag force errors at least for 2 nd , 3 rd and 4 th order Lagrange extrapolations. To conclude on the distance between Eulerian interpolation points and particle surface, as soon as δ ≥ ∆x, the error on drag force calculation is less than 2% for 3 rd and 4 th order Lagrange extrapolation in Stokes and Navier-Stokes regimes. In this case, a 3 rd order Lagrange extrapolation is considered, as being a good compromise between implementation complexity, calculation time and accuracy.

As illustrated in figure 9, the Lagrangian mesh refinement does not affect a lot the drag force computation as soon as the Lagrangian particle surface element size l is of the order of ∆x. For Re = 100, the error is around 2% on the drag force calculation.

In the rest of the present work, the 3 rd order Lagrange extrapolation will always be used together with l = ∆x.

Result on the drag coefficient

The transition area between axisymmetric flow and non-axisymmetric vortex shedding regime being around Re = 300, for a uniform flow past a fixed isolated sphere, we have conducted several simulations for Reynolds numbers up to 290 (Re = 0.1, 1, 10, 20, 40, 60, 80, 100, 150, 200, 250, 290). In each case, the simulations are stopped when the steady state of the flow is reached. The drag coefficient obtained at the final step of each simulation is compared to the correlations [START_REF] Beetstra | Drag force of intermediate Reynolds number flow past mono-and bidisperse arrays of spheres[END_REF] and [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF]. Figure 10 shows a very good agreement of the numerically calculated drag coefficient compared to the correlation of Schiller & Naumann [START_REF] Schiller | Über die grundlegenden Berechnungen bei der Schwerkraftaufbereitung[END_REF]. This test case provides an interesting validation of the force calculation for a single particle.

Pressure coefficient

The analysis of local pressure profiles can give insight into the pressure behavior depending on the Reynolds number. The local pressure coefficients are defined by equation ( 22) in a spherical coordinate system (see figure 11). They are considered along the polar angle θ as follows:

C p (θ) = p(θ) 1 2 ρU 2 ∞ ( 22
)
The pressure coefficient distribution according to θ is compared for Re = 1, 10, 100 to some available body fitted simulations results conducted by Magnaudet [START_REF] Magnaudet | Accelerated flows past a rigid sphere or a spherical bubble. Part 1. Steady straining flow[END_REF], Dennis & Walker [START_REF] Dennis | Calculation of the steady flow past a sphere at low and moderate reynolds numbers[END_REF], LeClair & al. [START_REF] Le Clair | A numerical study of the drag on a sphere at low and intermediate reynolds numbers[END_REF] and Massol [START_REF] Massol | Simulations numériques d'écoulements travers des réseaux fixes de spheres monodisperses et bidisperses, pour des nombres de Reynolds modérés[END_REF]. The It can be observed that a very good agreement is found between our force calculations and reference results. We can observe that the local pressure profile is symmetric with respect to θ = 90 for Re = 1 as shown in figure 12 a) and become increasingly asymmetric while increasing the Reynolds number in figures 12 b) and c). This feature of the flow is clearly highlighted for Re = 10 where the local pressure profile is no more symmetric with a negative pressure region that took place after θ = 60. For Re = 100, the same conclusion holds where the negative pressure region appears at θ = 50.

Forces in fixed arrangements of spheres

In the previous section, we have exhibited various successfull validations of the force calculation for uniform flows past a fixed isolated sphere. We can now confidently complexify the simulated particulate motion by investigating a flow past fixed particle assemblies which corresponds to a gas-solid flow with high Stokes number. Two way are possible to set up a uniform flow past fixed packed particles:

• by imposing a constant pressure gradient in the domain. In this configuration, the mean fluid velocity

u f = V (1 -C)udV V (1 -C)dV
evolves to reach a steady state corresponding to the imposed pressure gradient. • by choosing a desired mean fluid velocity u f = u d , and so a desired Reynolds number. A source term

F m = -ρ u f -u d ∆t + Np i F D i V f is in- serted in
Re S = (1 -α d )| u f |d ν ( 23
)
The relation between the solid volume fraction α d and the number of particle N p is:

N p = 6 π L D 3 α d (24)
As for the uniform flow past an isolated sphere, the grid resolution is fixed by imposing 5 Cartesian cells in the boundary layer as ∆x = d 5 √

Re S .

The mean non-dimensional drag force for all the particles is then defined as:

F = | F D | 3πµd(1 -α d )| u f | ( 25
)
with

F D = 1 N p Np i F D i , F D i
being the drag force computed over the i th particle.

Monodispersed arrangements of spheres

Face-Centred Cubic periodic arrangement of spheres

A Face-Centered Cubic (FCC) array is a cube where three spheres are placed on the faces centers, and one sphere is located on the vertices with periodic boundary conditions, as illustrated in the figure 13.

This configurations was widely studied [START_REF] Massol | Simulations numériques d'écoulements travers des réseaux fixes de spheres monodisperses et bidisperses, pour des nombres de Reynolds modérés[END_REF] [17] [START_REF] Tenneti | Drag law for monodisperse gas-solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres[END_REF] to understand and separate the wake effects, observed when spheres are aligned along the flow direction, i.e. streamwise interactions, from the blocking effects, where the spheres are aligned along the direction perpendicular to the flow direction, i.e. lateral interactions. Two regimes, i.e. attached and separated flows, govern the uniform flow past a FCC array of spheres, depending on the Reynolds number and the solid volume fraction. They are considered below.

Regime 1: attached flows

For low Reynolds numbers, the boundary layer remains attached to the particles. Unlike the isolated sphere configurations, the boundary layer detachment occurs for higher Reynolds numbers when FCC arrays are considered, due to the blocking effect of the surrounding particles. In fact, the presence of lateral spheres speeds-up the flow between the spheres and blocks the detachment of the boundary layer. Figure 14 illustrates this regime for Re = 50, α d = 0.15. A very good agreement is observed for the pressure distribution over the spheres in both plans β = 0 and β = 45 as shown in figure 14(a) and 14(b), compared to the results of Massol et al. [START_REF] Massol | Simulations numériques d'écoulements travers des réseaux fixes de spheres monodisperses et bidisperses, pour des nombres de Reynolds modérés[END_REF] obtained with body fitted meshes. The same agreement is observed for the axial friction coefficient as illustrated in figure With increasing Reynolds numbers, a separation of the boundary layer occurs in the downstream hemisphere of the particles. As noticed for an isolated sphere, the recirculating zone length increases with the Reynolds number. The flow confinement due to the presence of the lateral particles is so significant that high speed jet flows are obtained between the particles. These jet flows impact the spheres onto the upstream hemisphere. As illustrated in figure 15(a), this mechanism induces the development of a fountain effect on the vertical plan β = 0, the direct consequence of which is that the maximum pressure is no longer at θ = 0 but at θ = 36. Note that the same result was found by Massol et al. [START_REF] Massol | Simulations numériques d'écoulements travers des réseaux fixes de spheres monodisperses et bidisperses, pour des nombres de Reynolds modérés[END_REF]. Figures 15(a) and 15(b) show again the good agreement of the pressure distribution on the vertical and lateral plans with the body fitted simulations of Massol et al. [START_REF] Massol | Simulations numériques d'écoulements travers des réseaux fixes de spheres monodisperses et bidisperses, pour des nombres de Reynolds modérés[END_REF], and the same conclusions is observed in the figures 15(c) and 15(d) for the axial friction coefficient distribution. It can be noticed that for this Reynolds number of 300, larger differences are observed between our force calculation and body fitted grid results, mostly in the vicinity of the pressure peaks in the β = 0 plans.

Drag force

The drag force F normalized by the isolated sphere drag (given by Schiller [START_REF] Schiller | Über die grundlegenden Berechnungen bei der Schwerkraftaufbereitung[END_REF])

F s = 1 + 0.15Re 0.687 , noted as F F s
, is compared to existing literature results:

• Tenneti et al. [START_REF] Tenneti | Drag law for monodisperse gas-solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres[END_REF]:

F T = F s (1 -α d ) 2 + F φ + F φ,Re F φ = 5.81φ (1 -α d ) 2 + 0.48 α 1 3 d (1 -α d ) 3 F φ,Re = (1 -α d )α 3 d Re 0.95 + 0.61α 3 d (1 -α d ) 2
• Gobin et al. [START_REF] Gobin | Fluid dynamic numerical simulation of a gas phase polymerization reactor[END_REF]:

F G = F W Y if α d ≤ 0.3 min(F W Y , F E ) otherwise F W Y = F s (1 -α d ) -3.7
Wen&Yu [START_REF] Wen | Mechanics of fluidization[END_REF]]

F E = 150 18 α d (1 -α d ) 2 + 7 4 1 18 Re (1 -α d ) 2 Ergun[12]
• Beetstra et al. [START_REF] Beetstra | Drag force of intermediate Reynolds number flow past mono-and bidisperse arrays of spheres[END_REF]:

F B = 10α d (1 -α d ) 2 + (1 -α d ) 2 (1 + 1.5 √ α d ) + 0.413 (1 -α d ) 2 Re 24 1 (1-α d ) + 3(1 -α d )α d + 8.4Re -0.343 1 + 10 3α d Re -0.5-2α d (26)
Comparisons with these results are presented in figure 16 at different Reynolds number Re = 10 figure 16(a), Re = 50 figure 16(b), Re = 100 figure 16(c) and Re = 300 figure 16(d). It can be observed that a nice match is found between our results and those of Massol et al. [START_REF] Massol | Simulations numériques d'écoulements travers des réseaux fixes de spheres monodisperses et bidisperses, pour des nombres de Reynolds modérés[END_REF] which was extracted from body fitted simulations. Reasonable agreement is also observed with Gobin et al. [START_REF] Gobin | Fluid dynamic numerical simulation of a gas phase polymerization reactor[END_REF], and Beetstra et al. [START_REF] Beetstra | Drag force of intermediate Reynolds number flow past mono-and bidisperse arrays of spheres[END_REF] even if larger differences are noticed for Re = 10 and 300. The correlation given by Tenneti et al. [START_REF] Tenneti | Drag law for monodisperse gas-solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres[END_REF] seems to differ from the other results for high solid volume fraction. Note that the works of Beetstra et al. [START_REF] Beetstra | Drag force of intermediate Reynolds number flow past mono-and bidisperse arrays of spheres[END_REF] and Tenneti et al. [START_REF] Tenneti | Drag law for monodisperse gas-solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres[END_REF] was conducted using an immersed boundary methods to carry out these simulations.

Random periodic arrangement of spheres

A new step is taken in the complexity of the problem by randomly distributing a set of spheres for a given solid volume fraction. An illustration is given in figure 17 for Re S = 50 and α d = 0.3. The locations of the spheres are set by means of a classical Monte Carlo procedure to distribute them in the domain until the desired compacity is reached, and a hard sphere collisions model is used to prevent spheres from overlapping [START_REF] Campbell | Computer simulation of granular shear flows[END_REF]. To simulate an infinite domain, periodic boundary conditions are specified in all directions. To be statistically converged, we performed for most cases five Multiple Independent Simulations (MIS), Tenneti et al. [START_REF] Tenneti | Drag law for monodisperse gas-solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres[END_REF] conducted the same number of MIS to obtain their correlation for the drag force. Consequently, the representative non-dimensional drag force numerically obtained is an average of those calculated over the different realizations (MIS). Simulations with various solid volume fraction α d = 0.1, 0.15, 0.2, 0.3 have been investigated, together with different superficial Reynolds numbers ranging from 20 to 200.

The results presented in figure 18 demonstrate again a good global agreement of our results with existing correlations of Tenneti et al. [START_REF] Tenneti | Drag law for monodisperse gas-solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres[END_REF], Gobin et al. [START_REF] Gobin | Fluid dynamic numerical simulation of a gas phase polymerization reactor[END_REF] and Beetstra et al. [START_REF] Beetstra | Drag force of intermediate Reynolds number flow past mono-and bidisperse arrays of spheres[END_REF] for each solid volume fractions studied here as shown in figures: 18(a) for α d = 0.1, 18(b) for α d = 0.15, 18(c) for α d = 0.2 and 18(d) for α d = 0.3, all function of the superficial Reynolds number. In all cases, for larger Re S , the correlation of Beetstra et al. [START_REF] Beetstra | Drag force of intermediate Reynolds number flow past mono-and bidisperse arrays of spheres[END_REF] is farther from other works.

The pressure profiles according to θ are given in figure 20. They have been obtained with our force calculation and our simulations. No comparisons to other works is possible as these profiles are not reported in the literature. It can be noticed that the pressure contribution on the drag force increase with On the contrary, it seems that the Reynolds number has a small effect on the distribution of the pressure over the particles, and here again for both plans, as illustrated in figure 19. This result has previously been obtained in the study of Tenneti et al. [START_REF] Tenneti | Drag law for monodisperse gas-solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres[END_REF].

As can suggest the correlation of Gobin et al. [START_REF] Gobin | Fluid dynamic numerical simulation of a gas phase polymerization reactor[END_REF], our results (see figure 21) show that the drag force in a random arrangement of spheres is having the same dependence on the Reynolds number as for an isolated sphere. On the contrary, different behaviours are noticed for Tenneti and Beetstra simulations.

Bidisperse arrangements of spheres

For a bidisperse arrangement of sphere, i.e. with two types of particle size, we studied a Face-Centred Cubic periodic arrangement for the two species of spheres: the larger particles are distributed in the same configuration as the one previously presented for monodisperse arrangements while the smaller Two additional dimensionless parameters, for each species, are now necessary to characterize the flow:

x i = α i α d , y i = d i d s (27)
where d i and α i are the particle diameter and the solid volume fraction of the specie i respectively. The Sauter mean diameter d s is given by:

d s = 2 i x i d i -1
The Sauter mean diameter is one of the most important characteristic dimension for the bidisperse particle arrangement that is often used in the literature [START_REF] Beetstra | Drag force of intermediate Reynolds number flow past mono-and bidisperse arrays of spheres[END_REF]. 1.

Our force calculation results are given in figure 23. They are compared to Beetstra et al. [START_REF] Beetstra | Drag force of intermediate Reynolds number flow past mono-and bidisperse arrays of spheres[END_REF] correlation that provides the drag force for the i th species of particle as follows:

(F B ) i = ((1 -α d )y i + α d y 2 i + 0.0064(1 -α d )y 3 i )F B
where F B is the Beetstra drag force for a mono-dispersed arrangement of spheres given in equation [START_REF] Bizid | Penalty methods for turbulent flows interacting with obstacles[END_REF]. The values of forces obtained with our method are also compared to Massol body fitted results [START_REF] Massol | Simulations numériques d'écoulements travers des réseaux fixes de spheres monodisperses et bidisperses, pour des nombres de Reynolds modérés[END_REF]. It can be observed that for all Reynolds and solid fractions, our drag forces are in better agreement with Massol simulations than with Beetstra laws extracted from their simulations. A possible explanations of the gap between Beetstra correlation and other simulations is perhaps due to the accuracy and method for calculating forces on the particles with their approach. The same kind of discrepancies have been reported for Beelstra results in random arrangements of spheres.

Conclusions and perspectives

A new method has been designed for estimating forces in finite-size particle simulations. It is based on Lagrange extrapolation and Taylor interpolation of third order. In our method, a Lagrangian mesh is used to describe the particle shape over time on the Eulerian Cartesian mesh devoted to the flow motion. One of the most original part of our approach is to upwind the interpolation support in the normal direction to the fluid/particle interface in order to use only fluid values to estimate forces. Various parameters of the method have been tested such as effect of interpolation and extrapolation orders or the size of the particle surface elements on the error observed on forces.

Based on a second order fictitious domain method using penalty techniques and augmented Lagrangian procedures for the incompressibility constraint, we have simulated various particulate flow motions ranging from incompressible flows around an isolated particle at various Reynolds numbers to flows across packed spheres under Faced-Centered Cubic, random and bi-disperse arrangements. In all configuration, the drag forces have been compared to reference results of the literature for various solid fractions and Reynolds numbers. The general conclusion that we have obtained is that our force calculation method fits always better to body fitted simulations of Massol et al. [START_REF] Massol | Simulations numériques d'écoulements travers des réseaux fixes de spheres monodisperses et bidisperses, pour des nombres de Reynolds modérés[END_REF], that can be considered as a reference, compared to other force estimates obtained with other IBM approaches. It has also been demonstrated that having at least 5 points in the boundary layer attached to a particle ensures to recover all the physics of the interaction between the carrier fluid and particles without using any force or interaction model. To our understanding, for obtaining a correct force calculation in finite-size particle simulations, the most important parameter is not the number of Eulerian cells along a particle diameter but the number of mesh cells belonging to the boundary layer surrounding the particle, whatever the particle diameter.

Future works and ongoing developments are devoted to extending our method to heat flux calculation and extraction of heat transfer coefficients in particulate flows. interesting way of possible improvement of the accuracy of our approach is the use of Aslam extensions [START_REF] Aslam | A partial differential equation approach to multidimensional extrapolation[END_REF] to replace the Taylor interpolation in our force calculation method. Other interesting issues are the extension of the present force calculation to more complex and realistic particle shapes such as spheroid or spherocylindrical particles.

Acknowledgements

This work was granted access to the HPC resources of CINES under the allocation A0032b06115 made by GENCI (Grand Equipement National de Calcul Intensif), and to the resources of CALMIP supercomputing center under the allocation 2017-P1529."

Appendix 1: Taylor Interpolation

Let us consider f : X ∈ R 3 → R a n-differentiable function at a given point E. The n th order Taylor Interpolation of f at point C l is:

f (C l ) = n-1 |α|=0 1 α! ∂ α f ∂x α (E)(C l -E) α + O(||C l -E|| n )
where α = (α 1 , α 2 , α 3 ) ∈ N 3 is the sum multi-index. As described in figure 24, E denotes the nearest point to the Lagrangian particle surface point C l . Point E belongs to the fluid domain and so is located on the fixed Eulerian mesh.

The four first order or Taylor interpolation are the following:

• First order:

f (C l ) = f (E) + O(||C l -E||) (28)
• Second order:

f (C l ) = f (E) + ∂f ∂x i (E)((C l ) i -E i ) + O(||C l -E|| 2 ) ( 29
)
𝑛 𝑙 Fluid Solid 𝐶 𝑙 𝑆 𝐸 Fig. 24: Details of Lagrange extrapolation points for the drag force computation at point C l of the surface Lagrangian mesh. E is the nearest fluid point to the particle/fluid interface. The crosses represent the discrete compact support of the fluid points required for a 3 rd order Lagrange extrapolation.

• Third order:

f (C l ) = f (E) + ∂f ∂x i (E)((C l ) i -E i ) + 1 2 ∂ 2 f ∂x 2 i (E)((C l ) i -E i ) 2 + ∂ 2 f ∂x i ∂x j (E)((C l ) i -E i )((C l ) j -E j ) + O(||C l -E|| 3 ) (30)
• Fourth order:

f (C l ) = f (E) + ∂f ∂x i (E)((C l ) i -E i ) + 1 2 ∂ 2 f ∂x 2 i (E)((C l ) i -E i ) 2 + ∂ 2 f ∂x i ∂x j (E)((C l ) i -E i )((C l ) j -E j ) + 1 6 ∂ 3 f ∂x 3 i (E)((C l ) i -E i ) 2 + i =j 1 2 ∂ 3 f ∂x 2 i ∂x j (E)((C l ) i -E i ) 2 ((C l ) j -E j ) + ∂ 3 f ∂x 1 ∂x 2 ∂x 3 (E)((C l ) 1 -E 1 )((C l ) 2 -E 2 )((C l ) 3 -E 3 ) + O(||C l -E|| 4 ) (31)
For each order of interpolation, we have to discretize the partial derivative used in it by considering the same order of discretization. To avoid using solid point in the approximation of these derivatives, we choose to use, depending on the sign of the outgoing normal vector n l , a forward (resp. a backward) difference scheme if (n l ) i > 0 (resp. (n l ) i < 0). The forward difference scheme for the derivatives used in ( 28), ( 29), ( 30), ( 31) is detailed below. The backward difference scheme is straightforwardly obtained.

Gradient approximation

∂f ∂x i = 3f i -4f i+1 + f i+2 -2∆x i + O(||∆x|| 2 ) ∂f ∂x i = 11f i -18f i+1 + 9f i+2 -2f i+3 -6∆x i + O(||∆x|| 3 ) ∂f ∂x i = 25f i -48f i+1 + 36f i+2 -16f i+3 + 3f i+4 -12∆x i + O(||∆x|| 4 )
Hessian approximation

∂ 2 f ∂x 2 i = 35f i -104f i+1 + 114f i+2 -56f i+3 + 11f i+4 12∆x 2 i + O(||∆x|| 3 )
∂ 2 f ∂x i ∂x j = 539f i,j -781f i+1,j + 297f i+2,j -55f i+3,j -781f i,j+1 + 1035f i+1,j+1 192∆x i ∆x j -303f i+2,j+1 + 49f i+3,j+1 + 297f i,j+2 -303f i+1,j+2 + 3f i+2,j+2

192∆x i ∆x j + 3f i+3,j+2 -55f i,j+3 + 49f i+1,j+3 + 3f i+2,j+3 + 3f i+3,j+3 192∆x i ∆x j + O(||∆x|| 3 )

∂ 2 f ∂x 2 i = 45f i -154f i+1 + 214f i+2 -156f i+3 + 61f i+4 -10f i+5 12∆x 2 i + O(||∆x|| 4 )
∂ 2 f ∂x i ∂x j = 117f i,j -73f i+1,j -83f i+2,j + 42f i+3,j -3f i+4,j -73f i,j+1 -243f i+1,j+1 60∆x i ∆x j 477f i+2,j+1 -173f i+3,j+1 + 12f i+4,j+1 -83f i,j+2 + 477f i+1,j+2 60∆x i ∆x j + -528f i+2,j+2 + 137f i+3,j+2 -3f i+4,j+2 + 42f i,j+3 -173f i+1,j+3 60∆x i ∆x j + 137f i+2,j+3 -3f i+3,j+3 -3f i+4,j+3 -3f i,j+4 + 12f i+1,j+4 -3f i+2,j+4 60∆x i ∆x j + -3f i+3,j+4 -3f i+4,j+4 60∆x i ∆x j + O(||∆x|| 4 )

Third order derivative approximation

∂ 3 f ∂x 3 i = 49f i -232f i+1 + 461f i+2 -496f i+3 + 307f i+4 -104f i+5 + 15f i+6 -8∆x 3 i + O(||∆x|| 4 ) ∂ 3 f ∂x 2
i ∂x j = 2775f i,j -8725f i+1,j + 10895f i+2,j -7105f i+3,j + 2550f i+4,j -480∆x 2 i ∆x j -390f i+5,j -4085f i,j+1 + 11931f i+1,j+1 -13489f i+2,j+1 -480∆x + -34281f i,j+4,k+1 -7506f i+1,j+4,k+1 -6f i+2,j+4,k+1 -6f i+3,j+4,k+1 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+4,j+4,k+1 -7581f i,j,k+2 + 3444f i+1,j,k+2 -4356f i+2,j,k+2 -9000∆x 1 ∆x 2 ∆x 3 + -606f i+3,j,k+2 + 144f i+4,j,k+2 -136281f i,j+1,k+2 -9000∆x 1 ∆x 2 ∆x 3 + -4506f i+1,j+1,k+2 + 3744f i+2,j+1,k+2 -6f i+3,j+1,k+2 -6f i+4,j+1,k+2 -9000∆x 1 ∆x 2 ∆x 3 + 348519f i,j+2,k+2 + 3744f i+1,j+2,k+2 -6f i+2,j+2,k+2 -6f i+3,j+2,k+2 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+4,j+2,k+2 -257481f i,j+3,k+2 -6f i+1,j+3,k+2 -6f i+2,j+3,k+2 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+3,j+3,k+2 -6f i+4,j+3,k+2 + 51294f i,j+4,k+2 -6f i+1,j+4,k+2 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+2,j+4,k+2 -6f i+3,j+4,k+2 -6f i+4,j+4,k+2 -6f i,j,k+3 -9000∆x 1 ∆x 2 ∆x 3 + -1206f i+1,j,k+3 -6f i+2,j,k+3 + 494f i+3,j,k+3 -6f i+4,j,k+3 -9000∆x 1 ∆x 2 ∆x 3 + 97519f i,j+1,k+3 -6f i+1,j+1,k+3 -6f i+2,j+1,k+3 -6f i+3,j+1,k+3 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+4,j+1,k+3 -212981f i,j+2,k+3 -6f i+1,j+2,k+3 -6f i+2,j+2,k+3 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+3,j+2,k+3 -6f i+4,j+2,k+3 + 136469f i,j+3,k+3 -6f i+1,j+3,k+3 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+2,j+3,k+3 -6f i+3,j+3,k+3 -6f i+4,j+3,k+3 -20181f i,j+4,k+3 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+1,j+4,k+3 -6f i+2,j+4,k+3 -6f i+3,j+4,k+3 -6f i+4,j+4,k+3 -9000∆x 1 ∆x 2 ∆x 3 + -6f i,j,k+3 + 294f i+1,j,k+3 -6f i+2,j,k+3 -6f i+3,j,k+3 -6f i+4,j,k+3 -9000∆x 1 ∆x 2 ∆x 3 + -19531f i,j+1,k+3 -6f i+1,j+1,k+3 -6f i+2,j+1,k+3 -6f i+3,j+1,k+3 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+4,j+1,k+3 + 38294f i,j+2,k+3 -6f i+1,j+2,k+3 -6f i+2,j+2,k+3 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+3,j+2,k+3 -6f i+4,j+2,k+3 -18931f i,j+3,k+3 -6f i+1,j+3,k+3 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+2,j+3,k+3 -6f i+3,j+3,k+3 -6f i+4,j+3,k+3 -6f i,j+4,k+3 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+1,j+4,k+3 -6f i+2,j+4,k+3 -6f i+3,j+4,k+3 -6f i+4,j+4,k+3 -9000∆x 1 ∆x 2 ∆x 3 + O(||∆x|| 4 )

Fig. 1 :

 1 Fig. 1: Details of a 2D discretization of the particle surface S used in drag force computation

Fig. 3 :Fig. 4 :

 34 Fig. 3: Relative error (%) of the first component of velocity for Stokes flow past cylinder (left) and Stokes flow past sphere (right) in the whole domain

Fig. 5 :

 5 Fig. 5: Drag force relative error (%) according to extrapotation distance δ (written in ∆x unit) for the Stokes flow past a cylinder (left) and the Stokes flow past a sphere (right) at different order: ( ) first, ( ) second, ( ) third, and ( ) fourth orders are plotted.

Fig. 6 :

 6 Fig. 6: Drag force relative error (%) for the Stokes flow past a cylinder (left) and the Stokes flow past a sphere (right) in term of the Lagrangian mesh refinement l: force computed using 3 rd Lagrange extrapolation and 2 nd order Taylor interpolation

Fig. 8 :Fig. 9 :

 89 Fig. 8: Drag force relative error (%) for the uniform flow past sphere at Re = 100. Different orders of Lagrange extrapolation are considered: ( ) first, ( ) second, ( ) third, and ( ) fourth order. The distance between the first Eulerian point used to extrapolate forces and the particle surface is δ.

Fig. 10 :

 10 Fig. 10: Drag coefficient for uniform flow past sphere at different Reynolds number: ( ) Schiller & Naumann [11], ( ) present work.

Fig. 11 :Fig. 12 :

 1112 Fig. 11: Spherical coordinate system around a particle. The flow direction is represented by the undisturbed velocity U ∞ .(a) (b)

Fig. 13 :

 13 Fig. 13: Streamlines of a steady flow along the x-axis of a Face-Centred Cubic array of spheres at Re = 50 and α d = 0.5. The particle shape is plotted in blue.

Fig. 14 :

 14 Fig. 14: Pressure coefficient [(a) β 0, (b) β = 45] and axial friction coefficient [(c) β = 0, (d) β = 45], for a uniform flow past a FCC at Re = 50, α d = 0.15 : ( ) Massol et al. [23], and ( ) present work

Fig. 15 :

 15 Fig. 15: Pressure coefficient [(a) = 0, (b) β = 45] and axial friction coefficient [(c) β = 0, (d) β = 45], for a uniform flow past a FCC at Re = 300, α d = 0.15: ( ) Massol et al. [23], and ( ) present work

Fig. 16 :

 16 Fig. 16: Drag force for a uniform flow past a FCC, normalized by Schiller & Naumann [11] drag force for a uniform flow past an isolated sphere. Results are presented as a function of the solid volume fraction α d and Reynolds number (a)Re = 10, (b)Re = 50,(c)Re = 100 ,(d)Re = 300 : ( ) Tenneti et al. [22], ( ) Gobin et al.. [14], ( ) Beetstra et al. [19], ( ) Randrianarivelo [24], ( ) Massol et al. [23], and ( ) present work.

Fig. 17 :

 17 Fig. 17: Random periodic arrangement of spheres, for Re S = 50 and α d = 0.3 -The particle shape is plotted in blue while black lines represent streamlines.

Fig. 18 :

 18 Fig. 18: Normalized drag force for a uniform flow past random packed spheres as function of the superficial Reynolds number Re S and solid volume fraction (a) α d = 0.1, (b) α d = 0.15, (c) α d = 0.2, (d) α d = 0.3: ( ) Tenneti et al. [22], ( ) Gobin et al. [14], ( ) Beetstra et al. [19] and ( ) present work. The error bars represent 95% confidence intervals in the estimation of the normalized drag force.

Fig. 19 :Fig. 20 :

 1920 Fig. 19: Pressure coefficient for a uniform flow past random packed spheres at Re S = 200 as function of α d (a)β = 0, (b)β = 45: ( ) α d = 0.1, ( ) α d = 0.2 and ( ) α d = 0.3 (a) (b)

Fig. 21 :

 21 Fig. 21: Drag force for a uniform flow past random packed spheres, normalized by Schiller & Naumann [11] drag force for a uniform flow past an isolated sphere, as function of the superficial Reynolds number Re S and solid volume fraction (a) α d = 0.1, (b) α d = 0.15, (c) α d = 0.2, (d) α d = 0.3: ( ) Tenneti et al. [22], ( ) Gobin et al. [14], ( ) Beetstra et al. [19] and ( ) present work. The error bars represent 95% confidence intervals in the estimation of the normalized drag force

Fig. 22 :

 22 Fig. 22: FCC/FCC bidisperse arrangement of spheres for Re = 50 and α d = 0.3 -the particle shape is plotted in blue and the black lines are the streamlines.

Fig. 23 :

 23 Fig. 23: Normalized drag force for a uniform flow past a FCC/FCC packed spheres as function of the solid volume fraction α d (a) : (Re = 50, y1), (b) : (Re = 50, y2), (c) : (Re = 100, y1), (d) : (Re = 100, y2), : ( ) Beetstra & al. [19] , ( ) Massol [23], ( ) present work.

  Centered Cubic arrangement of spheres for both mono and bi-dispersed flows, and random assemblies of spheres only for mono-dispersed flows.

	The particulate Reynolds number Re = Centered Cubic arrangement of spheres following the work of Massol [23]. | u f |d is used studying the Face-ν
	On the other hand, in the random assemblies of spheres, another Reynolds
	number based on the superficial velocity (1 -α d )| u f | is considered to take into account the solid loading in the characterization of the flows [22][19]. It
	is given by:

the momentum conservation equations. It is adjusted until the desired Reynolds number is reached. This is the method used in this work. Its main advantage is to be able to simulate a prescribed Reynolds number without a trial and error procedure unlike to what is require in the pressure gradient technique.

Two ways to distribute the particles are studied in this work: a Face-

Table 1 :

 1 values of the parameters characterizing the bi-disperse flows simulated in this work. Reynolds number Re = 50 and 100. The solid volume fractions for the larger particles are α 1 = 0.15 and 0.3. The corresponding values of x i , y i and d s are reported in Table

	x 1		x 2	y 1	y 2	ds
	0.88889 0.11111 1.11111 0.55556 0.0009m
	We have performed simulations of FCC/FCC bidispersed arrangement of
	spheres, for which	d 1 d 2	= 2, at

  + f i+5,j+2 -385f i,j+3 + 711f i+1,j+3 -264f i+2,j+3 -64f i+3,j+3 -480∆x 2 i ∆x j + f i+4,j+3 + f i+5,j+3 + 10f i,j+4 -14f i+1,j+4 + f i+2,j+4 + f i+3,j+4 -480∆x 2 i ∆x j + f i+4,j+4 + f i+5,j+4 + 10f i,j+5 -14f i+1,j+5 + f i+2,j+5 + f i+3,j+5 ∂ 3 f ∂x 1 ∂x 2 ∂x 3 = 4069f i,j,k -13656f i+1,j,k -6f i+2,j,k -6f i+3,j,k -6f i+4,j,k -9000∆x 1 ∆x 2 ∆x 3 12194f i,j+1,k -23256f i+1,j+1,k + 10194f i+2,j+1,k + 444f i+3,j+1,k -9000∆x 1 ∆x 2 ∆x 3 + -306f i+4,j+1,k -16681f i,j+2,k + 77094f i+1,j+2,k -19206f i+2,j+2,k -9000∆x 1 ∆x 2 ∆x 3 + -456f i+3,j+2,k + 294f i+4,j+2,k -6f i,j+3,k -52956f i+1,j+3,k -9000∆x 1 ∆x 2 ∆x 3 + 11244f i+2,j+3,k + 244f i+3,j+3,k -6f i+4,j+3,k -6f i,j+4,k -9000∆x 1 ∆x 2 ∆x 3 + 13794f i+1,j+4,k -3006f i+2,j+4,k -6f i+3,j+4,k -6f i+4,j+4,k -9000∆x 1 ∆x 2 ∆x 3 + -6f i,j,k+1 + 18294f i+1,j,k+1 + 744f i+2,j,k+1 -6f i+3,j,k+1 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+4,j,k+1 + 59869f i,j+1,k+1 -6f i+1,j+1,k+1 -6f i+2,j+1,k+1 -9000∆x 1 ∆x 2 ∆x 3 + -6f i+3,j+1,k+1 -6f i+4,j+1,k+1 -177281f i,j+2,k+1 -40506f i+1,j+2,k+1 -9000∆x 1 ∆x 2 ∆x 3 + -756f i+2,j+2,k+1 -6f i+3,j+2,k+1 -6f i+4,j+2,k+1 + 153019f i,j+3,k+1-9000∆x 1 ∆x 2 ∆x 3 + 26994f i+1,j+3,k+1 + 1494f i+2,j+3,k+1 -6f i+3,j+3,k+1 -6f i+4,j+3,k+1 -9000∆x 1 ∆x 2 ∆x 3

			2 i ∆x j
	+ +	7911f i+3,j+1 -2654f i+4,j+1 + 386f i+5,j+1 + 1675f i,j+2 i ∆x j -480∆x 2 -3889f i+1,j+2 + 2856f i+2,j+2 -744f i+3,j+2 + 101f i+4,j+2 -480∆x 2 i ∆x j
	+	f i+4,j+5 + f i+5,j+5 -480∆x 2 i ∆x j	i ∆x j -480∆x 2 + O(||∆x|| 4 )


Valentin Gillot
email: valentin.gillot@irisa.fr

Frédéric Bimbot
email: frederic.bimbot@irisa.fr

Polytopic reconfiguration: a graph-based scheme for the multiscale transformation of music segments and its perceptual assessment

HAL is

INTRODUCTION

Music is usually considered as a sequential process, where sounds, group of sounds and motifs are occurring chronologically, following the natural unfolding of time. Under this approach, music is considered as a flow of information, the organization of which is essentially governed by the relations that develop between adjacent elements. This property of sequentiality is indeed central to many models

Copyright: c 2018 Valentin GILLOT of music description and generation, including musicological models such as Schenkerian analysis [START_REF] Schenker | Der freie Satz[END_REF], generative theories based on tree-like structures such as GTTM [START_REF] Lerdahl | A generative theory of tonal music[END_REF], or computational automata such as Markov chains [START_REF] Pearce | The construction and evaluation of statistical models of melodic structure in music perception and composition[END_REF].

From a complementary perspective, it is also commonly accepted that the organization of music is widely based on alternations of similarity and novelty which create patterns and systems of expectation and surprise that ultimately contribute to the structure of the musical discourse from a cognitive point of view [START_REF] Narmour | Music expectation by cognitive rulemapping[END_REF] [5] [START_REF] Bimbot | System & Contrast : A Polymorphous Model of the Inner Organization of Structural Segments within Music Pieces[END_REF]. These can be represented as a Polytopic Graph of Latent Relations [START_REF] Louboutin | Polytopic Graph of Latent Relations: A Multiscale Structure Model for Music Segments[END_REF] where each node of the graph represents a low-scale musical segment and vertices correspond to their relation within the expectation systems.

The aim of this work is to explore the polytopic model in terms of its relevance to account for the inner structure of musical segments, by assessing its potential use for the structural recomposition of music. In fact, as presented in section 2, the polytopic representation of a musical segment enables a large range of transformations by applying various permutations to its nodes, thus generating multiple reconfigurations of its musical content, with the same elements in a different order. Specific permutations, called Primer Preserving Permutations (PPP), are of particular interest, as they preserve systems of analogical implications between metrically homologous elements within the segment [START_REF]Description of Chord Progressions by Minimal Transport Graphs Using the System & Contrast Model[END_REF]. The central hypothesis of the present work is that the musical consistency of PPP-transformed segments will therefore be less affected than it would be by an "ordinary" (i.e. randomly generated) permutation.

In section 2, we describe the implementation of the polytopic reconfiguration process and we elaborate on the organizational properties of Primer Preserving Permutations as well as their potential impact on the inner structure of musical segments.

Then, in section 3, we assess the relevance of the reconfiguration scheme (and its underlying hypotheses): we report on a perceptual test where subjects are asked to rate musical properties of MIDI segments, some of them have been reconfigured with PPPs while others were transformed by Randomly Generated Permutations (RGPs) designed so as to possess a comparable number of discontinuities.

BACKGROUND AND FORMALISM

Musical objects and time-scales

Before detailing the background and formalism aspects of this work, it is important to specify what musical objects are under consideration in our study. Following Snyder's "levels of musical experience" [START_REF] Snyder | Music and memory: An introduction[END_REF], one may identify three ranges of time-scales in music: an event-fusion level, up to 1/32 th second (corresponding to early processing), an intermediate "melodic and rhythmic grouping" level, between 1/16 th and 8 sec. (governed by short-term memory) and a form level, 16 sec. or above (resorting to long-term memory).

In this work, we focus on sectional units (i.e. the first time-scale of the form level) and we study their inner organization as the result of the relations existing between successive constituents at the intermediate level, around 1 sec (i.e. typically half-bars). We put a particular focus on melody, with the consequence that our low-scale objects are small groups of a few notes with a common harmonic background.

Polytopic representations of implication systems

The basic concept of a sequential implication system is that the observation of a given event by a subject with no particular preconception will trigger (to that subject) the expectation that the next event is likely to be similar to the first one. This is considered as even more so, in the case of repetition, as is expressed by Narmour's implication principles for the analysis of basic melodic structures [START_REF] Narmour | The analysis and cognition of melodic complexity: The implication-realization model[END_REF]:

A + A -→ A A + B -→ C
thus meaning that the repetition of two similar patterns A induces the expectation of a third similar pattern A, whereas two different patterns A and B trigger the expectation of something different from both: C. Narmour's principles can be understood as a cognitive model based on Gestalt Theory, as mentionned in [START_REF] Morgan | Statistical learning and gestalt-like principles predict melodic expectations[END_REF].

This sequential implication principle can be extended by now considering a system of elements viewed as forming the base of an analogical induction process [START_REF] Bimbot | Semiotic structure labeling of music pieces: concepts, methods and annotation conventions[END_REF]. For instance, consider the matrices below, and imagine what could be the missing elements in terms of "logical" implications:

A B 2 4 ∩ ↓ E . 8 . ∪ .
As developed in [START_REF] Bimbot | System & Contrast : A Polymorphous Model of the Inner Organization of Structural Segments within Music Pieces[END_REF], such square systems trigger an expectation process which spans over 2 time-scales simultaneously and which writes, in a formalism compatible with Narmour's:

A + f (A) + g(A) -→ g(f (A))
and in particular:

A + B + A -→ B A + A + B -→ B
These two prototypical cases of square implication systems are indeed the basis of two very frequent structural patterns, ABAC and AABC, which can be understood as the denial, in 4 th position, of the implication system triggered by the first 3 elements. But quite a number of other such redundant patterns can be embedded in this framework.

Further generalizing this principle by encompassing more time-scales, as in [START_REF] Deruty | Methodological and musicological investigation of the System & Contrast model for musical form description[END_REF], leads to cubic (3-scale), tesseractic (4-scale) and more generally n-cubic systems where n time-scales are considered simultaneously (see Figure 1). When mapped with a sequence of musical events, these polytopic graphs can be used to represent analogical relations between musical objects in metrically homologous positions at different time-scales [START_REF] Louboutin | Polytopic Graph of Latent Relations: A Multiscale Structure Model for Music Segments[END_REF].

Structural reconfiguration by node permutations

There are many ways to define a permutation, but the most adequate one in the context of this work is to view it as a bijective function ϕ (i.e. a one-to-one correspondence) between the N time indexes of an original sequence x 0 ...x N -1 and those of the permuted sequence x ϕ(0) ... x ϕ(N -1) .

Applying a permutation to the time indexes creates a certain degree of disruption in the original flow of the sequence. In this paper, we characterize this fact by two properties: D, the number of discontinuities created by the permutation and E, the (log 2 of the) maximum timeinterval (or excursion) between 2 consecutive elements in the permuted sequence.

D = #{t | [ϕ(t) -ϕ(t -1)] = 1} 0<t<N E = log 2 max 0<t<N |ϕ(t) -ϕ(t -1)|
Among all possible permutations ϕ, a particular subset of them is focused on in this work. They are referred to as PPPs (for Primer Preserving Permutations) and were introduced by Louboutin et al. [START_REF]Description of Chord Progressions by Minimal Transport Graphs Using the System & Contrast Model[END_REF] as the set of permutations which preserve the systemic relations of the musical elements x t in the sequence, by just interchanging the time-scales at which they develop. These permutations are therefore particularly interesting to investigate on the relevance of the implications systems (and their preservation) in the consistency of a musical segment.

In the case of a sequence of N =16 elements, the polytope is a 4-cube, also called a tesseract (see Figure 2), and it can be viewed as being composed of 4 homologous lowerscale systems of 4 elements (formed by 4 parallel faces of the tesseract), themselves linked by an upper-scale system formed by the first element of each lower-scale system (on a face which is perpendicular to the 4 lower-scale ones) (see Figure 3). Mapping these systems of faces by a permutation, while preserving time order within faces, yields 5 possibilities (plus the identity), which correspond to the PPPs defined in Table 1. Note that each PPP maps 0 to 0 and N -1 to N -1. By construction, PPPs preserve implication systems of 4 elements in the sense that sets of 4 parallel faces in the tesseract are mapped on another set of parallel faces. By doing so, analogical systems are preserved, and only the scales at which they develop are modified: for instance,

[x 0 x 1 x 2 x 3 ] in the original is "transferred" to [x 0 x 4 x 1 x 5 ]
with PPP1 and [x 0 x 4 x 8 x 12 ] with PPP4. Moreover, the last element of a system never occurs before all the elements of all the systems it belongs to, have themselves occurred, therefore preserving a principle of causality in all the analogical implications.

For comparison purposes, Table 2 illustrates 5 Randomly Generated Permutations (RGPs), which also start with 0, end with N -1, and have the same profile as their corresponding PPPs in terms of discontinuity and excursion. It can be easily checked that, unless incidentally, analogical systems are not preserved by RGPs.

Ultimately, a given permutation is applied to a musical se-quence of events by reading and copying (in the permuted order) the original musical material after having segmented into proper low-scale elements. Figure 4 illustrates this process on a 8-bar melodic line (original at the top), processed as a sequence of N =16 half-bar (2-beat) elements, resulting in reconfigured melodic lines, the organization of which is either driven by PPPs or by RGPs.

Handling musical surface time-shifts

As a prerequisite to the application of a permutation to a sequence of musical objects, it is of primary importance to define precisely the location and boundaries of the objects that will be subjected to the transformation. In this work, elementary objects are melodic fragments However, for reconfiguration purposes, these "objects" may not be optimal if they are blindly synchronized on downbeats (and sub-downbeats). In a number of cases, the musical surface of the "proper" low-scale melodic elements should be slightly shifted in order to correspond to an optimal grouping of the notes and therefore a better result after permutation. This situation is illustrated in Figure 5 where the first melodic element should be considered as starting 1.5 beat before the first downbeat of the section.

In designing the perceptual test described in the next section, this optimal time-shift has been adjusted manually on a case-by-case basis. However, we present in the appendix, a preliminary algorithm that has been developed to estimate the time-shift automatically.

PERCEPTUAL EXPERIMENTS IN STRUCTURAL RECONFIGURATION

In this section, we present the details of the methodology which has been used for the evaluation of polytopic reconfigurations. The general principle has been to design a perceptual test where subjects are presented with musical excerpts which are either PPP-transformed or RGPtransformed, and where they are asked to rate (blindly) these excerpts both in terms of consistency and attractive-ness. We then study the difference of ratings across subjects between the two types of transformations.

Test data

The experiments are performed on the RWC POP corpus, a subset of the Real World Computing (RWC) music database created for scientific purposes by Goto et al. [START_REF] Goto | RWC Music Database: Popular, Classical and Jazz Music Databases[END_REF]. This dataset has been designed for researchers. It is available for a small fee and without restrictions.

The RWC POP set is made of 100 songs in WAV and MIDI formats. According to the authors, they were composed partly in the style of the 80s American hits and partly in the 90s Japanese style. This corpus is very commonly used in various task in Music Processing and Information Retrieval.

Using the RWC songs in our tests has an additional advantage: there is virtually no risk that subjects participating to the test already know the songs from which the excerpts are stemming. This is a favorable situation as the familiarity with a particular song could affect the rating of the transformed excerpts in an uncontrolled way.

For our experiments, a subset of 24 songs has been selected for which the first instance of chorus is exactly 8bar long. The reconfiguration process is applied to the MIDI version of the chorus. The elements of interest are the melodic lines, the accompaniment (harmony and bass) and the drums. We consider that these elements provide an acceptable rendering of the structural information of the chorus and they all are transformed synchronously by the permutations. The average duration of a chorus is 19.5 seconds.

As presented in section 2, these transformations are performed on the original chorus after having segmented it into 16 elements (with a possible shift of the musical surface of the melody, as explained in section 2.4). Then, the 10 permutations (5 PPPs and 5 RGPs) are applied to the 16 elements, resulting in a re-ordering of the original choruses, i.e. 10 new musical contents (per chorus), which we will call "variants" of the original.

To summarize, for each of the 24 choruses, the following 11 variants are considered:

• the original version of the chorus

• 5 PPP-transformed variants (PPP 1 to PPP 5 ) • 5 RGP-transformed variants (RGP 1 to RGP 5 )
The perceptual test is therefore conducted on 24 × 11 = 264 distinct excerpts.

Protocol

The perceptual experiments are performed with 66 subjects, selected among colleagues and relatives of the authors 1 . About 25% of the subjects are female and 75% male, aged between 22 and 55 years old.

In order to avoid the bias that could arise if a subject gets familiar with a particular chorus, each subject is presented with one and only one variant of each of the 24 choruses, in a random order. The total number of tests in our experiments is therefore equal to 66 × 24 = 1584, i.e. 144 tests per variant. Moreover, the 11 variants are distributed as evenly as possible for each subject, i.e. typically 2 (or occasionally 3) per subject.

After a brief general presentation of the test (during which the subjects are not provided with any detail on the scientific motivation of the test), they are asked to answer, for each excerpt, two questions formulated as follows:

1. "Please give your opinion as of the consistency of the musical construction of this passage".

"Please give your opinion as of the degree of musical attractiveness of this passage"

The subjects must respond by placing a cursor on a 5level scale corresponding to 5 different levels of evaluation (from very bad consistency to very good consistency and from very unattractive to very attractive). They also have the possibility to position the cursor between two graduations, thus offering them 9 levels of rating altogether.

Asking two different questions aims at leading subjects to decorrelate as much as possible their rating of the acceptability of the musical construction (question 1), from their own personal taste (which they can express in their response to question 2).

Scoring

For each subject, we thus collect 24 consistency ratings and 24 attractiveness ratings.

Let us denote as r, the generic variable which designates such ratings (ignoring for the moment whether it is consistency or attractiveness, for the sake of notations' simplicity).

The 24 ratings r jk produced by subject j (1 ≤ j ≤ 66, 1 ≤ k ≤ 24) are distributed into the 11 variants: O for original, P i for PPP i (1 ≤ i ≤ 5) and P * i for RGP i (1 ≤ i ≤ 5). These individual ratings are averaged separately, thus leading to 11 mean ratings (one per type of variant):

µ 0j = 1 m0j
O r jk for the original choruses [START_REF] Schenker | Der freie Satz[END_REF]

µ ij = 1 mij
Pi r jk for the 5 PPP variants

(2)

µ * ij = 1 m * ij P * i r jk for the 5 RGP variants (3)
Note that, in practice, m 0j , m ij and m * ij are equal to 2 or 3 and their sum is equal to 24.

To evaluate the relative degradation created by the various reconfigurations of the music material, we can measure:

δ ij = µ 0j -µ ij (4) δ * ij = µ 0j -µ * ij (5)
namely the signed difference, for a given subject j, between his/her ratings of the reconfigured variants and his/her ratings of the originals.

For each subject and variant, we can also compute:

π ij = µ ij -µ * ij = δ ij -δ * ij ( 6
)
which measures the (positive or negative) preference of PPP i -transformed excerpts over RGP i -transformed ones, for speaker j (where PPP i and RGP i are comparable in terms of discontinuities D and excursion E).

Ultimately, we focus on 3 scores:

π j = 1 5 5 i=1 π ij (7)
δ i = 1 66 66 j=1 δ ij (8)
δ * i = 1 66 66 j=1 δ * ij ( 9
)
where π j provides a global preference score (positive or negative) on PPPs over RGPs for each subject j, while δ i (resp. δ * i ) provides a degradation score for each of the PPP-(resp. RGP-) reconfigured variants, averaged over all subjects.

These two quantities are calculated separately for consistency ratings and for attractiveness ratings.

Results

Figure 6 and 7 depict the distributions of consistency and attractiveness preference scores (π j ), ranked from lowest to highest, over the panel of 66 subjects. In both cases, the distributions are visibly shifted towards positive values, especially for consistency ratings. Indeed, for consistency, 57.5 subjects have a positive preference for PPPs versus 8.5 having a negative preference 2 (i.e. 87.1 % vs 12.9 %). The average score for positive judgments amounts to +0.89 whereas it rests at -0.30 only, for negative judgments. For attractiveness, the proportions are still clearly in favor of PPPs, but not as contrasted (74.2 % vs 25.8 %), with average scores of +0.74 and -0.55 respectively.

Figure 8 and 9 represent the relative degradation scores (δ i and δ * i ) for the 10 variants, ranked in increasing order of degradation. From the ranking of the permutations, it is noticeable that the degradation score is primarily correlated with the number of discontinuities of the permutations, as the two permutations with lower D show a smaller value of δ as opposed to the three others. Globally, these results indicate that PPP-based reconfigurations tend to be less disruptive than RGPs against the 2 counting as 0.5 in each category, the subject who scores exactly 0. perceived consistency of musical segments, thus supporting the hypothesis that the preservation of expectation systems contributes to their inner structure.

CONCLUSIONS

Systems of relations between elements play a key role in structural constructions in many cognitive dimensions. The experiments reported in this paper constitute an initial investigation on the role of analogical systems (modeled by polytopic graphs) in the perception of segmental structure in pop music.

These results indicate a very noticeable trend towards the hypothesis that the preservation of the analogical implication systems in musical segments impacts positively the perception of their inner structure and hence, their acceptability. However, more extensive tests on a larger population and with more varied musical excerpts are certainly needed to consolidate these first results, and refine their precise scope.

In parallel, we intend to develop the proposed reconfiguration scheme into a music creation tool which could be used by composers to experiment new possibilities in employing musical material with multiple purposes, results and effects. In this context, the polytopic representation is an additional advantage, as it is bound to make manipulation interfaces more intuitive and easy to use.

Ultimately, the development of an automatic segmentation process is a complementary goal towards a fully operational concept.

Appendix Melodic segmentation with a compressibility criterion

In order to perform a relevant permutation-based transformation to a sequence of musical objects, it is crucial to define precisely the location and boundaries of the objects on which the transformation is applied.

As mentioned in section 2.4, it is quite common that the structure of the melodic line is time-shifted with respect to the phase of the metrical pulsation. As a consequence, it is important to take into account this time-shift when defining the low-level melodic objects. In this section, we present a preliminary method to automatically estimate this timeshift.

For doing so, we estimate a segmentation of the melodic surface by optimizing the grouping of the notes using a compressibility criterion.

Meter, segmentation and phase-shift

For simplicity and consistency with the rest of our work, we consider musical segments which are 8-bar long, and 4|4 meter, i.e L = 32 beats, and we assume that we want to segment the excerpt in m elementary objects. In this work m = 16.

Let S be the sequence of instants where the strong beats fall, according to the meter and the bar:

S = (t 0 , t 1 , ..., t p , ..., t m-1 ) (10)
Let us now denote as Σ another segmentation that is not necessarily synchronized with downbeats. Σ = (τ 0 , τ 1 , ..., τ p , ..., τ m-1 )

In the general case we can write:

τ p = t p + θ p (12)
where θ p is the phase-shift of the musical surface of segment p in relation to the strong beat location t p .

In compliance with the work reported in the body of this article, we further assume that θ p does not depend on p:

τ p = t p + θ (13)

Properties of the phase-shift

When θ = 0, the melodic segments are considered to be starting synchronously with the downbeat. If θ > 0, this means that the melodic objects display some phase-delay with the beat, whereas if θ < 0, they start in anticipation (as in the case of an anacrusis).

Let α be the "size" of one elementary object in terms of the number of beats, i.e.:

α = L m (14)
Here α = 2.

Quite naturally, we assume that the phase-shift cannot exceed the size of the elementary objects, i.e.:

-α ≤ θ ≤ α (15)
Fig. 10, shows the segmentation of a 8-bar melody, represented in MIDI format. The first 3 notes of this section are located with a phase-shift of -1,5 (beat) with respect to the metric boundaries: the thick lines correspond to the downbeat instants and the thin ones correspond to the segmentation obtained with θ = -1.5 (beat).

A compressibility criterion for segmentation

In order to estimate an optimal value of θ, we define a compressibility criterion which scores how redundant is the sequence of elements resulting from a given segmentation.

This approach is inspired from Kolmogorov's information theory [START_REF] Kolmogorov | Three approaches to the quantitative definition ofinformation[END_REF] and relates, in its spirit, to recent work exploring the potential of approaches based on complexity and compression for modeling music contents (for instance [16] [17]).

Let now Σ be a segmentation as in equation 11 and σ p the elementary object of Σ within the time interval [τ p , τ p+1 [. We first define the parallel run-length between two melodic elements as the relative duration d ∈ [0, 1] during which they remain similar to each other (i.e. identical, up to a translation). We denote as b the translation between the two similar fragments and we assign to b a binary value of 0 (if no translation) or 1 (when a translation is observed). We then define an elementary cross-compressibility score function c(σ p , σ q ) between two melodic objects σ p and σ q , by combining d and b: c(σ p , σ q ) = d(σ p , σ q ) -λ b(σ p , σ q ) (16)

For the time being, λ is tuned empirically.

We then define a compressibility score function for each elementary object σ p within Σ, which depends on the previous objects in the segmentation, for instance:

z(σ p ) = 1 p h<p c(σ h , σ p ) (17)
Finally, we compute the entire segmentation compressibility score Z(Σ) as:

Z(Σ) = 1 m p z(σ p ) ( 18
)
Because θ is the same across segments, all possible segmentations Σ can be tested exhaustively, and the optimal segmentation is chosen as the one with the highest compressibility score.

Current performance level

In its current version, the results provided by this algorithm have been evaluated over the 24 choruses used in the perceptual experiments, by comparing the automatic estimations of the melodic phase-shifts with those determined by a manual annotation: in 10 cases (out of 24), the automatic algorithm provided the same result as the expert.

This approach is currently being improved but the current results encourage us to further investigate the method, so as to evaluate its full potential.

Figure 1 .

 1 Figure 1. Graph representation of multi-scale implication systems. From left to right: segment, square, cube, tesseract.

Figure 2 .

 2 Figure 2. Example of the polytopic representation of a time sequence in the case of a tesseract.

Figure 3 .

 3 Figure 3. Representation of the 6 systems of parallel faces on the tesseract and their corresponding PPP in reference to the original.

Figure 4 .

 4 Figure 4. Example of the 11 types of reorganized variants on a 8-bar melodic line.

Figure 5 .

 5 Figure 5. Example of a melodic line starting 1.5 beat before the first downbeat of the section.

Figure 6 .

 6 Figure 6. Musical consistency: relative PPP-preference πj of the test subjects, ranked in increasing order

Figure 7 .

 7 Figure 7. Musical attractiveness: relative PPP-preference πj of the test subjects, ranked in increasing order

Figure 8 .Figure 9 .

 89 Figure 8. Musical consistency: degradation scores δi and δ * i for the different permutation variants (PPP and RGP)

Figure 10 .

 10 Figure 10. Example of a melodic line segmentation starting 1.5 beat before the first downbeat of the section (θ = -1.5 beat).

Fig. 11

 11 depicts the behavior of Z for a particular example, with θ ranging from -2 to 2 by steps of 0.25 beat. Score values are normalized to 1 for the maximum value.

Figure 11 .

 11 Figure 11. Behavior of the compressibility score Z for a particular 8-bar melody, with θ ranging from -2 to 2 beats.

Table 1 .

 1 Extensive definition of the 6 PPPs and their properties (discontinuity and excursion -see text).

	PPP

Table 2 .

 2 Examples of RGPs and their properties (discontinuity and excursion -see text).

For obvious reasons, the authors themselves were not included in the panel of subjects.