Adrien Scheuer
email: adrien.scheuer@uclouvain.be.

Amine Ammar

Emmanuelle Abisset-Chavanne

Elías Cueto

Francisco Chinesta

Roland Keunings

Suresh G Advani

Elias Cueto

Data-driven upscaling of orientation kinematics in suspensions of rigid fibres

Keywords: Fibre suspensions, data-driven upscaling, closure approximations

come

Introduction

In processes involving fibre suspensions ( e.g. composite manufacturing, papermaking, biological and pharmaceutical applications, food-processing and cosmetics industries, etc.), predicting the evolution of particle orientation is critical since the rheology of the material and its final properties depend on the microstructure. Classically, three modelling scales can be distinguished: the microscopic, mesoscopic and macroscopic scales.

At the microscopic scale, the orientation of a single particle is identified by a unit vector p aligned with the particle axis. In the case of a Newtonian suspending fluid, the evolution of the rod orientation is governed by Jeffery's equation [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF]]

ṗJ = ∇v • p -(∇v : (p ⊗ p))p, (1)
whith ∇v the unperturbed fluid velocity gradient. These kinematics, derived under the assumptions of a Stokes flow, lay the foundation for nearly all models used today. Extending Jeffery's theory to account for various internal or external effects, including Brownian effects [START_REF] Chinesta | From single-scale to two-scales kinetic theory descriptions of rods suspensions[END_REF]], bending phenomena [Abisset-Chavanne, Férec, Ausias et al. (2015)], particle inertia [START_REF] Scheuer | Modelling the effect of particle inertia on the orientation kinematics of fibres and spheroids immersed in a simple shear flow[END_REF]], electrical forces [START_REF] Perez | On the multi-scale description of electrical conducting suspensions involving perfectly dispersed rods[END_REF]], wall effects [START_REF] Perez | A multi-scale description of orientation in simple shear flows of confined rod suspensions[END_REF]; Scheuer, Abisset-Chavanne, Chinesta et al. (2016)] is readily achievable using a dumbbell representation of a rod [START_REF] Bird | Dynamics of Polymeric Liquids[END_REF]; [START_REF] Binetruy | Flows in Polymers, Reinforced Polymers and Composites: A Multiscale Approach[END_REF]]. Despite the richness of the possible descriptions at the microscopic scale, the computational effort to efficiently track millions of particles (as in scenarios of industrial interest) is in general unaffordable. Coarser descriptors are thus called for. At the mesoscopic scale, the information regarding the orientation state of a population of particles is contained in a scalar probability density function (pdf) ψ(x, t, p), that provides the fraction of particles with a given conformation p at any position x and time t. Solving the associated Fokker-Planck equation, governing the time evolution of the pdf, is however a challenge for traditional numerical methods, due to the inherent high-dimensionnality of the problem. Particle methods have long been used to conduct simulations at that scale [Öttinger (1996)], and very few works addressed the continuous Fokker-Planck equation [START_REF] Lozinski | A fast solver for Fokker-Planck equation applied to viscoelastic flows calculations: 2D FENE model[END_REF]; [START_REF] Chauvière | Simulation of dilute polymer solutions using a Fokker-Planck equation[END_REF]]. Notable progress were made recently [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF]] with the introduction of the Proper Generalized Decomposition [START_REF] Ammar | A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids[END_REF], 2007)], able to address high-dimensional PDEs.

At the macroscopic scale, the pdf is substituted by its first moments, provinding a crude, yet concise description of the orientation state in the material. In the case of fibres, due to the symmetry of the pdf, odd-order moments vanish. The so-called second and fourth order orientation tensors, introduced by Advani et al. [START_REF] Advani | The use of tensors to describe and predict fiber orientation in short fiber composites[END_REF]], read respectively

a(x, t) = S (p ⊗ p)ψ(x, t, p)dp (2)
and

A(x, t) = S (p ⊗ p ⊗ p ⊗ p)ψ(x, t, p)dp, (3)
where the integration is performed on the unit sphere S. These orientation tensors exhibit particular properties. a and A are completely symmetric, that is

a ij = a ji (4)
and

A ijkl = A ijlk = A kijl = A lijk = A klij , etc. (5)
Moreover, due to the normalisation condition of ψ, the trace of a equals 1.

The second-order orientation tensor has an intuitive physical interpretation. A high value of the first (resp. second and third) diagonal component of a indicates that the particles tend to orient along this x-(resp. yand z-) direction. If all the diagonal components are 1 3 , the 2 orientation tensor suggests three-dimensional random orientation state, but triaxial or any other orientation state that gives this average is also possible. This is an example of the inherent ambiguity of crude macroscopic descriptors. When two diagonal components are equal to 1 the tensor suggests two-dimensional random or planar biaxial orientation state. Finally, a unit diagonal component indicates full alignment in that direction.

The time evolution of the second-order orientation tensor is readily obtained using Jeffery's kinematics Eq. ( 1)

ȧ = ∇v • a + a • (∇v) T -2A : ∇v. (6)
However, this expression involves the fourth-order orientation tensor A. Unfortunately, the time derivative of the fourth-orientation tensor, using the same rationale, involves the sixth-orientation tensor and so on. Thus, a closure approximation is required. Much research has focused on developing accurate and stable closure approximations, indicating that the problem is far from being solved. We propose in the sequel an overview of the closures proposed in the literature; an in-depth discussion of the subject can be found in Jack et al. [Jack andSmith (2004, 2007)].

• Simple closures: linear (LIN) [START_REF] Hand | A theory of anisotropic fluids[END_REF]] (exact in the case of isotropic orientations), quadratic (QUAD) [START_REF] Hinch | Constitutive equations in suspension mechanics. Part 2: Approximate forms for a suspension[END_REF]] (exact for aligned fibres) and hybrid (HYBR) [START_REF] Advani | The use of tensors to describe and predict fiber orientation in short fiber composites[END_REF]] (combining the two previous ones);

• Composite closures: attempting to approximate directly the second-order tensor A : ∇v [START_REF] Hinch | Constitutive equations in suspension mechanics. Part 2: Approximate forms for a suspension[END_REF]];

• Orthotropic closures: attempting to express A in the principal axis of a [START_REF] Cintra | Orthotropic closure approximations for flow-induced fiber orientation[END_REF]];

• Natural closures: natural (NAT) [START_REF] Dupret | Modeling the flow of fiber suspensions in narrow gaps[END_REF]] and IBOF [START_REF] Chung | Invariant-based optimal fitting closure approximation for the numerical prediction of flow-induced fiber orientation[END_REF]] are fitted closures based on the most general expansions of A in terms of a and ∇v;

• Neural-newtork-based closures: NNET [START_REF] Jack | Neural network-based closure for modeling short-fiber suspensions[END_REF]] and NNORT [START_REF] Qadir | Modeling fibre orientation in short fibre suspensions using the neural network-based orthotropic closure[END_REF]];

• Closures for the sixth-order orientation tensor: such as LIN 6 , QUAD 6 , HYBR 6 [START_REF] Advani | Closure approximations for three-dimensional structure tensors[END_REF]], or even invariant-based fitted closures INV 6 and IBF 6 [Jack andSmith (2005, 2006)].

The macroscopic scale offers a simple and crude description of the microstructure. Simulations at that scale are thus much cheaper, explaining why this description is preferred in industrial applications. The pdf is substituted by some of its moments, sacrificing the level of detail and the involved physics in favour of computational efficiency. Closure approximations remain however an issue.

In this work, we propose a methodology aimed at providing data-driven macroscopic simulations of orientation kinematics that are cheap and closure-free. The approach consists of an offline step, the construction of a database of scenarios obtained from accurate microscopic simulations, and an online step, the data-driven macroscopic simulation itself. Data-driven approaches have surged over the last decade as a new paradigm in simulationbased engineering. Unprecedented possibilities were introduced in Dynamic Data Driven Application Systems (DDDAS) [START_REF] Darema | Dynamic data driven applications systems: A new paradigm for application simulations and measurements[END_REF]; [START_REF] Michopoulos | Dynamic-data-driven real-time computational mechanics environment[END_REF]], entailing the ability to incorporate additional data into an executing application, improving modelling methods and prediction capabilities. [START_REF] Kirchdoerfer | Data-driven computational mechanics[END_REF]] developed a strategy to carry out mechanical calculations directly from experimental material data (and pertinent constraints and conservation laws), thus bypassing the empirical material modelling step of conventional approaches. Closely related, the works of Peherstorfer et al. [Peherstorfer andWillcox (2015, 2016)] focus on constructing reduced-order models directly from data, inferring the full-order operators without the need to construct them explicitly, or even to have a direct knowledge about the governing models. The paper is structured as follows. Section 2 explains the main idea of our data-driven upscaling approach. In Section 3, this methodology is first illustrated in the well-known case of dilute fibre suspensions, where it can be compared against macroscopic closurebased models. Its relevance is then shown in the case of confined fibre suspensions, for which closures proved to be inadequate [START_REF] Perez | A multi-scale description of orientation in simple shear flows of confined rod suspensions[END_REF]; Scheuer, Abisset-Chavanne, Chinesta et al. (2016)]. Finally, we apply this framework in a more complex case involving semi-concentrated suspensions of electrically-charged rods, for which no reliable macroscopic model is available. We draw in Section 4 the main conclusions of this work.

2 Data-driven upscaling of orientation kinematics

The main idea behind our data-driven approach is the following: Since the physics at the microscopic scale can be modelled reasonably enough, we can conduct expensive accurate offline direct numerical simulations at that scale and extract the corresponding macroscopic descriptors in order to build a database of scenarios. During the online stage, the macroscopic descriptors can then be updated quickly by combining adequately the items from the database instead of relying on a sometimes imprecise macroscopic model (usually involving closure approximations). This methodology is depicted schematically in Fig. 1. Specifically, the two stages are as follows:

• Offline stage: construction of the data-base. A large amount of microscopic simulations involving populations of N fibres are run, exploring a wide range of initial orientation configurations. Moreover, the different scenarios may also include variations in the suspension parameters, such as the flow velocity gradient, the fibre volume fraction (influencing the inter-particle interactions), the confinement state (see below), the applied electric field in the case of charged particles (see below); these parameters are collectively referred to as α. For each population of fibres, the macroscopic descriptors, (a, a ˙, α), are computed from the microscopic ones, (p i , p ˙i, α), i = 1, . . . , N , in order to build a database of scenarios. An additional step, not explored in this paper, is to reconstruct a map a ˙= f (a, α) from interpolations of the items in the database [Gonzálea, Chinesta and Cueto (2018)].

(p i , ṗi , ↵) (a, ȧ, ↵) ȧ = f (a, ↵) Offline stage Online stage or f (a, ↵) a(t, ↵) a(t + t, ↵) = a(t, ↵) + t ȧ ȧ t t + t
Figure 1: Data-driven approach to fibre orientation kinematics Our aim is to present a general framework for the upscaling of orientation kinematics in suspensions of rigid fibres. More advanced and sophisticated techniques can be used to perform efficiently both offline and online stages. On the one hand, the offline construction of the database can be recast in the general context of Design of Experiments (DoE) methodology [START_REF] Eriksson | Design of Experiments[END_REF]], to help decide which scenarios of interest to consider. Latin hypercube sampling [START_REF] Mckay | Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF]] is a popular sampling technique to explore parameter spaces in experimental designs. On the other hand, during the online data-driven simulation, the comprehensive lookup in the database could be substituted either by the evaluation of an interpolation map f (as discussed in the perspectives) or alternatively using classical regression methods over the database items, or using machine-learning techniques and neural-networks trained over the database entries.

• Online stage: macroscopic data-driven simulation. At each time step, we identify in the database the closest items to the current orientation state (and suspension parameters) a(t, α) and combine them to obtain the instantaneous evolution kinematics. In the case where the mapping f was built, this evolution is readily obtained using f . And so on for the next time steps.

In this section, we propose an illustration of the methodology that has just been presented, in the case of suspensions of rods. In particular, we explain how to practically construct the database and the metrics use to measure the distance between orientation tensors. We first discuss the classical unconfined dilute case, for which there is a long history of macroscopic models, allowing us to assess the performance of our approach. We will then discuss the relevance of the method for confined suspensions, for which traditional macroscopic models fail. Finally, we will discuss the case of semi-concentrated suspensions of electrically-charged fibres, using microscopic direct numerical simulations inspired by molecular dynamics.

3.1 Unconfined dilute suspensions of rods immersed in a simple shear flow

In the case of dilute suspensions of rods immersed in a Newtonian fluid, the kinematics of each particle follow Jeffery's equation, Eq. ( 1). In this section, we consider that the suspension undergoes a simple shear flow, whose velocity field is given by v = γz 0 0 T , with γ = 1 s -1 . In this flow, it is known that the fibres simply align in the flow field (x-direction).

Database construction

The database is built by following the evolution of populations of N = 5000 particles.

In order to cover a wide range of initial configurations, a physically-reasonable choice is to consider initial orientations as "Gaussian" distributions. Since fibre orientations can be depicted as points on the unit sphere, we consider Von Mises-Fisher distribution of mean m (unit orientation vector) and variance s. Fig. 2a depicts an example of such a distribution. These orientation states were sampled from the Von Mises-Fisher distributions using the normal-tangent decomposition property of the distribution on the sphere, see [START_REF] Chen | Generate random samples from von mises-fisher and watson distributions[END_REF]] for additional details and the Matlab implementation. In this study, we construct two databases, the first with (n m , n s ) = (10, 15) (150 initial configurations), and the second, more comprehensive with (n m , n s ) = (20, 12) (240 initial configurations). The n m individual mean values are uniformly distributed over the sphere and the n s variances range from 0.05 (fibres nearly aligned) to 1.75 (fibres nearly uniformly distributed all over the sphere). For each population, we run the flow simulation during 30 seconds, and compute each 10 -1 second the macroscopic descriptors a and ȧ from the individual p i and ṗi

(i = 1, . . . , N ) as a = 1 N N i=1 p i ⊗ p i (7) and ȧ = 1 N N i=1 ṗi ⊗ p i + p i ⊗ ṗi . (8)

Data-driven simulation

During the online stage, we identify in the database the K items a db k (k = 1, . . . , K) closest to the current orientation tensor a(t). To do so, we use the euclidean distance applied to the vectorized forms of the orientation tensors, composed of its independent components:

vec(a) = a 11 a 22 a 12 a 13 a 23 T . (9)
Then, we compute a weighted average of the K corresponding kinematics ȧdb k to derive the instantaneous evolution ȧ that can be applied to have the orientation tensor at the next time step, that is a(t + ∆t) = a(t) + ∆t ȧ.

(10)

The reconstruction weights are obtained by solving the minimization problem

min wk vec a(t) - K k=1 w k vec a db k 2 (11) such that K k=1 w k = 1 (12) w k ≥ 0. ( 13
)
Remark: The choice of the definition of distance is definitely a de licate question. An ideal choice would be to have access to the "geodesic" distance on the manifold described by the trajectories of the second-order orientation tensors, but such a distance is far from obvious. In this work, we choose to stick with the Euclidean distance (as described above), that provided satisfactory results, as long as there are enough samples on the manifold. Fig. 3 shows two examples of simulations. In each case, only the diagonal components of the orientation tensors are depicted: the solid colour lines correspond to the discrete orientation tensor (computed for validation purposes); the discontinuous colour lines correspond to the data-driven orientation tensor (here K = 5) and the discontinuous grey lines correspond to the closure-based macroscopic models (for the QUAD, HYBR and IBOF closures). The solid colour lines correspond to the discrete approach (computed for validation purposes); the discontinuous colour lines correspond to the data-driven approach and the discontinuous grey lines correspond to the closure-based macroscopic models (for the QUAD, HYBR and IBOF closures)

In these examples, we can see that the data-driven simulations perform quite well, as does the macroscopic model using the fitted IBOF closure. The QUAD and HYBR closures tend however to accelerate the orientation transients.

Regarding the computational costs, closure-based models run two order of magnitude faster (0.03 s) than microscopic simulations (30 s using N = 5000 fibres). The data-driven ap-proach lies in-between, requiring 3s, but there is room for improvements since we consider here, as a proof-of-concept, a naive implementation (using extensive searches in the database to identify the neighbouring points for example).

Performance assessment

In order to assess properly the performance of our approach, we compare the predictions of the macroscopic data-driven simulations and closure-based macroscopic models with microscopic simulations. Specifically, we average the L2 relative error computed on the first diagonal component a 11 of the second-order orientation tensor over n c random initial distributions composed of N fibres. The average L2 relative error is defined as

Ē = 1 n c nc c=1 E c , (14)
with

E = a macro 11 (t) -a micro 11 (t) 2 dt a micro 11 (t) 2 dt , (15)
where a macro (t) is computed either using the data-driven approach or a closure-based model. As before, we consider the QUAD, HYBR and IBOF closure approximations and a micro (t) is obtained form the expensive discrete microscopic simulations. The error estimator here is based on the first diagonal component of the orientation tensor only, since in the case of a shear flow, the suspending fibres tend to align in the flow field, and thus a 11 is an approximate indicator of the alignment state as it approaches unity.

The results of this comparative study (n c = 500) are shown in Fig. 4. We observe that the data-driven approach shows improved performance compared to conventional closurebased models (QUAD or HYBR) but state-of-the-art fitted closures (IBOF) still provides the lowest upscaling error. As expected, using the most comprehensive database (with (n m , n s ) = (20, 12)) improves the accuracy of the data-driven method.

Confined dilute suspensions of rods immersed in a simple shear flow

We now move to confined suspensions of rods, that is suspensions flowing in gaps narrower than the fibre length. The gap walls now prevent the particle from rotating freely and some trajectory, passing outside the flow domain are thus forbidden. We have shown in previous work [START_REF] Perez | A multi-scale description of orientation in simple shear flows of confined rod suspensions[END_REF]; Scheuer, Abisset-Chavanne, Chinesta et al. (2016Chinesta et al. ( , 2018))] that in that case, the fibre kinematics can be written as Jeffery's equation augmented with an additional term that prevents the fibre from leaving the flow domain, ṗ = ṗJ + ṗC ,

where ṗC = -

( ṗJ •n) (1-(p•n) 2 ) (n -(p • n)p)
, with n a unit vector normal to the gap wall (the contact force is assumed to be orthogonal to the wall, friction being neglected). Equipped with these new fibre kinematics, we can follow the same rationale as in the unconfined case. We consider a configuration where the confinement is strong: the ratio between the gap height H and the fibre length L is set t o H = 0 .2. The suspension undergoes the same simple shear flow as before.

Database construction

The construction of the database is similar as in the previous case, except that the initial orientation states are now given by distributions that are Gaussian in the azimuthal direction and uniform across the narrow gap height. An example is depicted in Fig. 2b. The mean vectors m are uniformly distributed on the equator, and the variance s ranges from 0.05 (concentrated) to 1.75 (all over the allowed domain).

Data-driven simulation The data-driven simulation proceeds exactly as described in the unconfined case. Fig. 5 shows two examples of simulations. In each case, only the diagonal components of the orientation tensors are depicted: the solid colour lines correspond to the discrete orientation tensor using the confined kinematics Eq. ( 16) (computed for validation purposes); the discontinuous colour lines to the discrete orientation tensor using Jeffery's kinematics Eq. ( 1) (computed to assess the impact of confinement); the discontinuous dashed colour lines correspond to the data-driven orientation tensor (here K=5) and the discontinuous grey lines correspond to the closure-based macroscopic models (for the QUAD, HYBR and IBOF closures). In both examples, we note that the closure-based macroscopic models completely fail to address confinement configurations, even when the impact of confinement on the kinematics itself is low (in situations where few fibres actually interact with the gap wall as in Fig. 5, bottom). In other words, and as concluded in our previous work [START_REF] Perez | A multi-scale description of orientation in simple shear flows of confined rod suspensions[END_REF]; Scheuer, Abisset-Chavanne, Chinesta et al. (2016)], the main challenge with traditional macroscopic models involving moments of the orientation pdf lies more with representation capabilities in highly confined conditions than with a suitable description of the induced orientation kinematics. On the other hand, the data-driven approach reproduces quite well the predictions provided by the expensive microscopic simulations.

Confined kinematics

Figure 6: Average L2 relative upscaling error for the macroscopic models either computed using the data-driven approach (DB) or a closure approximation (QUAD, HYBR, IBOF):

(left) comparison against the unconfined kinematics Eq. ( 1); (right) comparison against the confined kinematics Eq. ( 16)

Performance assessment We use the same method as before to assess the performances of the method (here n c = 100). The results are depicted in Fig. 6. This comparative study confirms the observations of the previous figure (Fig. 5). As shown in Fig. 6 (right), traditional closure-based models (even using a robust fitted closure) tend to mispredict the orientation kinematics by more than 15%, whereas the data-driven approach still concedes only 5% of relative error. For the sake of completeness, Fig. 6(left), computes the relative error with respect to the (hypothetical) unconfined kinematics, to support our claim that closure approximations are inadequate for initial confined configurations (independently of the kinematics itself).

Dilute suspensions of rods immersed in a complex flow

If we consider that the particles are immersed in a complex flow (instead of a simple shear flow), the same rationale can be applied. Indeed, in the dilute regime, the fibre kinematics are governed by Jeffery's kinematics Eq. ( 1), which shows a linear dependency with the velocity gradient. Thus, databases can be built for elementary flows (for example: shear flow in the x-, yand z-directions, uniaxial elongation in the xand y-directions and rotation flow around the x-, yand z-directions), and during the online stage, the local velocity gradient is decomposed in its elementary contributions, and the outcomes of the different databases are weighted accordingly.

Semi-concentrated suspensions of electrically-charged rods

In the remainder of this section, we address the kinematics of electrically-charged rods (dipoles) immersed in a Newtonian fluid and subject to an external electric field E. A multi-scale modelling of such suspensions (in the dilute case) was already proposed in [START_REF] Perez | On the multi-scale description of electrical conducting suspensions involving perfectly dispersed rods[END_REF]]. A microscopic model governing the evolution ṗ of a single rod is obtained from a micromechanical derivation using a dumbbell representation of the particle and reads ṗ = ṗJ + ṗE , where ṗE depends on the external electric field E and the charge q of the rod dipole. The proposed macroscopic model is however tainted with non-reliable closure approximations that make it impractical to use. We could of course apply the same rationale again, building on top of this modified Jeffery equation. However, we want to consider here a semi-concentrated suspension, that is we want to account for the effects of fibre-fibre interactions as well. Moreover, we want to emphasize that the data-driven methodology proposed in this work is general and does not depend on the technique used to conduct the microscopic simulations. Therefore, we use microscopic direct numerical simulations inspired by molecular dynamics (MD). This fine scale s imulation technique i s based o n the following a ssumptions: (i) each rod consists of a set of connected particles; (ii) inter-particle interactions are described from appropriate potentials, in particular, the Lennard-Jones potential V LJ and two other potentials, V E and V B , used to describe respectively the rod elongation and bending; (iii) the rods are subject to inertial, hydrodynamic (drag) and electrical forces. A description of the inner workings of this molecular dynamics simulation is out of the scope of this paper but the details can be found in Perez et al. [START_REF] Perez | On the multi-scale description of micro-structured fluids composed of aggregating rods[END_REF]]. Specifically, the microscopic MD simulations follow the evolution of N electrically-charged rods interacting with each other in a periodic representative volume element. Fig. 7a shows the initial isotropic configuration of the particles. As before, we consider here a simple shear flow, whose velocity field is given by v = γz 0 0 T , with γ = 1 s -1 . The electric field points upwards in the z-direction and the charge q on each rod extremity is set to q = 1 C. In the absence of an electric field, the fibres tend to align in the flow field, as illustrated in Fig. 7b, that shows the final orientation state of the fibres when E = 0 NC -1 . Conversely, when the electric field is strong, the fibres cannot align in the flow and the final orientation is along an inclined axis (the inclination depends on the intensity of the electric field), as illustrated in Fig. 7b, that shows the final orientation state of the fibres when E = 50 NC -1 .

Database construction

In this illustration, we only vary the number of particles N in the suspension (that is the concentration of the suspension and thus the potential number of inter-particles interactions) and the intensity of the external electric field E. Adding a variation of the shear rate γ or the initial orientation state (as in two previous illustrations) is a straightforward extension. Therefore, the databases are built by following the evolution of populations of N = 100, 200, . . . 800 particles subjected to an electric field E of intensity ranging from 0 to 60 NC -1 . The shear rate is fixed, γ = 1 s -1 and the isotropic orientation state is always chosen as initial configuration. In each case, we run the MD simulation during 10 seconds, and compute each 10 -1 second the macroscopic descriptors a and ȧ from the individual p i and ṗi (i = 1, . . . , N ). Due to the stochastic nature of the MD simulations, the simulations are run ten times and the results are averaged over the ten realizations. Data-driven simulation The data-driven simulation is a bit different from what was presented in the two previous illustrations, since we now have two parameters that influence the kinematics of the suspension: The number of particles in the system N (that influences the amount of fibre-fibre interactions) and the intensity of the external electric field E. During the online stage, we intend to carry a simulation characterized at each time t by the current number of fibres in the system (N t ) and the current value of the electric field intensity (E t ). Among the databases at our disposal, we then identify the ones that best match the value of the current parameters (for example if (N t , E t ) = (225, 35) we keep the four databases built for N = 200 and 300 and E = 30 and 40 NC -1 ) and compute the weights needed for a bilinear interpolation of these results. For each one of these, we proceed as described before, looking within the individual database to find the K closest orientation tensors to the current orientation tensor and combining them adequately. Finally, these individual results are then combined using the bilinear weights computed just before. These manipulations might appear a bit tedious but are fairly easy and they actually provide a flexible way to handle for example time-varying electric fields fr om th e static databases. Proceeding in this way allows us to actually interpolate among the parameter space, even though the parameters are not of the same order of magnitude. Indeed, interpolating directly on the data triplet (N, E, a) would not have provided meaningful results, at least with the Euclidean distance, due to disparity of the quantities at stake. Fig. 8 shows three examples of simulations, in the case of weak, medium and strong external electric field. In each case, only the diagonal components of the orientation tensors are depicted: the solid colour lines correspond to the discrete orientation tensor obtained from MD simulations (computed for validation purposes) and the discontinuous colour lines correspond to the data-driven orientation tensor. As described previously, in the case of a nearly zero electric field ( Fig. 8, top), the fibres tend to align in th e flow field (xdirection) and thus the first diagonal component of the orientation tensor is dominant. On the contrary, when the electric field is strong ( Fig. 8, bottom), the particles are mostly aligned in the z-direction and thus the third diagonal component of a is important. When the electric field is of medium intensity ( Fig. 8, middle), the fibres tend to align in an intermediate orientation and the first and third components of a are in balance. In the three examples, the data-driven approach was in excellent agreement with the fine-scale simulations.

Performance assessment Again, we use the same method as before to assess the performances of the method. The number of random configurations (value of N and E ) is n c = 100. In this case, there is however no macroscopic model to compare with. We found that the data-driven method concedes only 5.9% of relative error with respect to fine-scale MD simulations. Regarding the computational costs, in this example, the data-driven approach runs in less than a second whereas MD simulations, inherently expensive, require from 30 to 500 seconds depending on the number of particles in the system.

Conclusion and perspectives

We presented a data-driven methodology aimed at providing efficient closure-free macroscopic simulations of the orientation of suspended rigid fibres, using a database of pre- The solid colour lines correspond to the discrete approach obtained from MD simulations (computed for validation purposes) and the discontinuous colour lines correspond to the data-driven approach. From top to bottom: weak, medium and strong external electric field E computed scenarios obtained from accurate direct computations at the microscopic scale. We show the relevance of this approach in the well-known case of dilute fibre suspensions, where it performs as well as state-of-the-art closure based models, but also for suspensions of confined or electrically charged fibres, for which conventional closure-based methods proved to be inadequate and reliable macroscopic models are simply not available. Therefore, this method appears as an appealing and "easy-to-set-up" technique in situations where closure-based models are unsatisfactory or have not been developed yet, including in situations where the physics at stake is complex (for example in the case of fibre-fibre interactions), provided that adequate microscopic simulation techniques are available.

In addition to the many situations where this methodology could be applied, many perspectives are envisioned for a data-driven approach in the context of fibre suspensions. First, even at the microscopic scale, Jeffery's equation could be replaced by some kinematics learned from experimental observations, especially in the case of non-Newtonian matrix suspensions for which there is no counterpart available (with the exception of Brunn's work [START_REF] Brunn | The slow motion of a rigid particle in a second-order fluid[END_REF]] for second-order fluid in the limit of low Weissenberg and the recent multi-scale modelling based on that model [START_REF] Borzacchiello | Orientation kinematics of short fibres in a second-order viscoelastic fluid[END_REF]). Second, a data-driven approach to predictions of the suspension rheology is to be explored. Another track, mentioned in the description of our data-driven approach but not explored in this paper, is the possibility of interpolating the items in the databases in order to build an approximation map that could be used directly during the online stage. The recent works on multi-dimensional interpolation techniques based on the Proper Generalized Decomposition (PGD), in particular the sparse PGD [Ibáñez, Abisset-Chavanne, Ammar et al. (2018a)] and the local PGD [Ibáñez, Abisset-Chavanne, Chinesta et al. (2018b)], open the way for interesting perspectives in that direction.

  Confined case: Gaussian in the azimuthal direction and uniform across the gap height

Figure 2 :

 2 Figure 2: Example of distributions of initial orientations used to build the database

Figure 3 :

 3 Figure3: Evolution of the diagonal components of the orientation tensor a for unconfined dilute suspensions of rods. The solid colour lines correspond to the discrete approach (computed for validation purposes); the discontinuous colour lines correspond to the data-driven approach and the discontinuous grey lines correspond to the closure-based macroscopic models (for the QUAD, HYBR and IBOF closures)

Figure 4 :

 4 Figure4: Average L2 relative upscaling error for the macroscopic models either computed using the data-driven approach (DB) or a closure approximation (QUAD, HYBR, IBOF)

Figure 5 :

 5 Figure5: Evolution of the diagonal components of the orientation tensor a for unconfined dilute suspensions of rods. The solid colour lines correspond to the discrete orientation tensor using the confined kinematics Eq. (16) (computed for validation purposes); the discontinuous dash-dotted colour lines to the discrete orientation tensor using Jeffery's kinematics Eq. (1) (computed to assess the impact of confinement); the discontinuous dashed colour lines correspond to the data-driven orientation tensor (here K=5) and the discontinuous grey lines correspond to the closure-based macroscopic models (for the QUAD, HYBR and IBOF closures)

  Final orientation state in the absence of an electric field (E = 0 NC -1 ): the fibres align with Final orientation state in the presence of a strong electric field (E = 50 NC -1 ): the fibres align along an intermediate inclined axis

Figure 7 :

 7 Figure 7: Periodic representative volume element containing N = 500 interacting electrically-charged fibres used in the molecular dynamics simulations. The suspension is subject to a simple shear flow and an electric field E is applied in the +z-direction

Figure 8 :

 8 Figure8: Evolution of the diagonal components of the orientation tensor a for semi-concentrated suspensions of electrically-charged rods. The solid colour lines correspond to the discrete approach obtained from MD simulations (computed for validation purposes) and the discontinuous colour lines correspond to the data-driven approach. From top to bottom: weak, medium and strong external electric field E

ICTEAM, Université catholique de Louvain, Av. Georges Lemaitre 4, Louvain-la-Neuve B-1348, Belgium.

Arts et Métiers ParisTech, Boulevard du Ronceray 2, BP 93525, 49035 Angers cedex 01, France.

ICI & ESI GROUP Chair, Ecole Centrale de Nantes, Rue de la Noe 1, F-44300 Nantes, France.

Aragon Institute of Engineering Research, Universidad de Zaragoza, Edificio Betancourt. Maria de Luna, s.n.,

Zaragoza, Spain. 5 PIMM & ESI GROUP Chair, ENSAM ParisTech, Boulevard de l'Hôpital 151, F-75013 Paris, France. 6 Deptartment of Mechanical Engineering, University of Delaware, Newark, DE 19716, UnitedStates.

Acknowledgements: A. Scheuer is a Research Fellow of the "Fonds de la Recherche Scientifique de Belgique" -F.R.S.-FNRS.


Pierre Puchaud

Georges Dumont

Nicolas Bideau

Charles Pontonnier

Knee Torque Generation Capacities Modelled With Physiological Torque-Angle-Velocity Relationships

Keywords: maximal joint torque, isokinetic dynamometer, torque-angle-velocity relationship, knee, muscle mechanics

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The torque generation capacities of a human body are fundamental factors governing human movement. They are often assessed for human performance, as well as in musculoskeletal modelling. In such models, scaling individual strength generation capacities is challenging but can provide physiologically meaningful data. Basically, models are fitted to isokinetic measurements of joint torques in different angle and angular velocity conditions. Assuming muscles are viscoelastic actuators, their entire architectures contribute to Joint Torque-Angle and Torque-Velocity Relationships (JTAR and JTVR respectively, and their coupling JTAVR) at the joint level. On one side, experimental observation at different scales (muscle sarcomere, muscle fibre and joint) resulted in various JTAR models available in the literature. On the other side, JVTR models are often modelled without obvious physiological consistency. Recently, a new JTVR model was proposed increasing physiological transparency of the elbow JTAVR [START_REF] Haering | Using Torque-Angle and Torque-Velocity Models to Characterize Elbow Mechanical Function: Modeling and Applied Aspects[END_REF]), but those results might be jointspecific. The current study aimed at evaluating five JTAR and two JTVR models on the knee flexion and extension of 27 participants at maximal intensity.

Materials and methods

Data collection

Twenty-seven healthy males (25.9 ± 4.5 years; 1.8 ± 0.1 m; 75 ± 11.6 kg) gave their consent for the study (CPP SOOM I, France no. 2018-A00484-51). Participants seated with their dominant leg alongside a Con-Trex MJ® isokinetic dynamometer (CMV AG, Dübendorf, Switzerland). The dynamometer axis was aligned with the knee epicondyles in flexion at 90°. Tight straps immobilized the thigh and the shank. Range of motion was adjusted to the subject. Goniometric measurements were used to calibrate angular values. After a warm-up, five maximal voluntary knee flexions and extensions hold for four seconds at angles equally spaced through the range of motion were recorded as isometric trials. Then, three repetitions of concentric isokinetic trials were recorded at 60°.s -1 , 120°.s -1 , and 180°.s -1 in flexion and extension. Participants were strongly encouraged to produce maximal torque values.

Joint Torque-Angle-Velocity Relationships

As previously described in [START_REF] Haering | Using Torque-Angle and Torque-Velocity Models to Characterize Elbow Mechanical Function: Modeling and Applied Aspects[END_REF], five JTAR were used as Quadratic, Cosine, Cubic, Sinusexponential, Normal models computing the isometric torque Γ for a given angle α, defined by:

Γ = Γ • JTAR(α, α , 𝑅𝑜𝑀) (1)
Isometric parameters of JTAR are the maximal isometric torque Γ , the optimal joint angle α and the maximal range of isometric force production 𝑅𝑜𝑀. Two JTVR were used as Anderson-based [START_REF] Anderson | Maximum Voluntary Joint Torque as a Function of Joint Angle and Angular Velocity: Model Development and Application to the Lower Limb[END_REF]) and Power-based [START_REF] Haering | Using Torque-Angle and Torque-Velocity Models to Characterize Elbow Mechanical Function: Modeling and Applied Aspects[END_REF]) to compute the concentric torque Γ for a given angular velocity ω. It was defined by:

Γ = Γ • JTVR(ω, 𝑷) (2)
For the Anderson-based model, concentric parameters 𝑷 were velocity at 75% of maximal isometric torque and ratio between velocities at 50% and 75% of maximal isometric torque. For the Power-based model, concentric parameters 𝑷 of were maximal concentric velocity and velocity at maximal power.

Parameters identifications

A least-square-curve-fitting method (trust-region algorithm, Matlab® Optimization Toolbox TM ) minimized the quadratic distance between modelled and experimental torques, optimizing successively isometric and concentric parameters α , 𝑅𝑜𝑀 and 𝑷 in a two-steps optimization.

Model Evaluation

Root mean square errors (RMSE) between modelled and experimental knee flexion/extension torques were computed for each JTAR models and for each combination of JTAR and JTVR models. Identified JTVR parameters of models were compared. A twoway ANOVA was performed between Andersonbased and Power-based RMSE for each JTAR to evaluate the effect of JTVR models on RMSE. Concerning the JTAR models, differences in fitting experimental data were observed in Table 1, Figure1. The lowest RMSE were displayed for the two symmetrical models (quadratic and cosine) for flexion and extension. The quadratic and cosine JTAR appeared as the most accurate and adaptable model on the large cohort, as for the elbow joint. Mean RMSE of cosine JTAR model were the lowest in extension. RMSE were higher because of higher experimental torques, compared to previous study on the elbow.

Results and discussion

Figure 1: Torque-Angle-models fitting of a typical participant for knee extension.

Concerning the JTVR models, differences in fitting were also observed: an example is provided

Conclusions

Five JTAR and two JTVR were compared when fitting dynamometric measurements of the knee. While a quadratic and cosine JTAR model fitted best the isometric data, the recently proposed Power-based JTVR increased physiological transparency for human performance without decreasing the data fitting.

*Corresponding author. Email: pierre.puchaud@irisa.fr

  Figure 2. knee flexion RMSE were lower than knee extension. The Anderson-based model displayed lower RMSE than the Power-based model in flexion and extension. Mean velocities at 75% of maximal isometric torque identified among participants were between 65.3 and 101.4 deg/s (cubic and sinus-exp) in flexion, and between 38.1 and 55.3 deg/s (normal and cubic) in extension.Mean ratios between velocities at 50% and 75% of maximal isometric torque were between 3.24 and 3.75 (sinus-exp and normal) in flexion, and between 3.9 and 4 (cosine and normal) in extension. The Power-based model displayed the highest RMSE in mean for flexion (not significant p>0.05) and extension (significant p<0.05). Mean maximal concentric velocities computed among participants were between 979.7 and 1455 deg/s (normal and sinusexp) in flexion, and between 702 and 853 deg/s (normal and cubic) in extension. Mean velocities at maximal power were between 251.9 and 448.7 deg/s (cubic and sinus-exp) in flexion, and between 175.6 and 213.5 deg/s (normal and cubic). The velocity at maximal power, also called optimal velocity is a widely used criterion of human performance[START_REF] Mcdaniel | Joint-specific power-pedaling rate relationships during maximal cycling[END_REF]. Thus, the Powerbased model provide an insightful model to link dynamometric measurements to human performance.

Figure 2 :

 2 Figure 2: Knee flexion/extension Torque-Angle-Velocity fitting of a typical subject with the quadratic JTAR and the Power-based JTVR. Black dots and lines correspond to isometric and isokinetic measures.

Table 1 :

 1 Knee flexion/extension mean RMSE between modelled and experimental maximal torques for the five JTAR and two JTVR models(Anderson-based and Power-based).

		Models Quadratic Cosine	Cubic Sinus-exp Normal
		JTAR	11.7±4.9 13.6±5.4 26.1±8.5 33.9±16.7 39.8±9.0
	Flex	Anderson 11.2±3.5 11.5±3.3 13.6±4.1 29.9±23.4 22.9±5.0
		Power	13.7±5.3 14.2±5.2 17.4±6.1 31.6±22.0 25.1±5.4
		JTAR	17.3±6.6 15.7±6.7 18.2±11.2 31.0±10.5 36.9±10.1
	Ext	Anderson 29.0±8.6 26.8±7.4 25.4±10.4 29.1±7.9 29.4±7.1
		Power	34.8±10.4 33.0±9.7 31.8±11.2 37.1±9.3 37.6±8.9


Fast marching and acoustic descriptors based method for fish proportion interpolation

Imen Karoui, Ronan Fablet and Jean-Marc Boucher Telecom Bretagne, UMR CNRS lab-STICC 3192, UEB, France Abstract-We propose a new method for the estimation of fish abundance from both acoustic data and some trawl hauls catches. In this work, we operate at a global level and we aim at estimating fish abundance from these images and not to identify the species of each school. We associate each trawl catch to the nearest acoustic image and we describe each image by a set of global statistical distributions estimated on it and related to each school parameters. Then, we use the fast marching algorithm, a front propagation based scheme to define a region of interest around each trawl associated image. The fast marching algorithm propagates each front initialized on the image associated to the trawl samples with a velocity proportional to the distance between the trawl image acoustic features and those of the image for which we want to estimate the fish abundance. Finally, the fish abundance of each image is estimated as a weighted sum of the abundances associated to each trawl. The weights are estimated from the propagation time given by the fast marching algorithm. The method is experimented on real and synthetic data. INTRODUCTION Fish abundance estimation is the cornerstone of ecosystem based fisheries. Trawl surveys are the most important fisheries data source used in stock assessment. The acoustic survey data, with its more resolved sampling structure, is used to establish a relationship between the acoustic and trawl data and to interpolate the fish abundance in zones where trawls are not available. Here, we aim at developing an automatic method for the estimation of the fish abundance. In the literature, we distinguish two main approaches:

1-School based approaches that extract school fish from the echograms and state the problem as the classification of these extracted fish schools [START_REF] Hammond | A general procedure for estimating the composition of fish school clusters using standard acoustic survey data[END_REF][START_REF] Haralambous | Artificial neural netorks as a tool for species identification of fish schools[END_REF]. 2-Global approaches that are stated at the echogram level and estimate the global fish abundance on a whole acoustic image (an echogram part) [START_REF] Petitgas | Sampling variance of species identification in fisheries-acoustic surveys based on automated procedures associating acoustic images and trawl hauls[END_REF][START_REF] Hammond | Bayesian estimation of fish school cluster composition applied to a Bering Sea acoustic survey[END_REF]. The most popular method is the algorithm AICASA (Acoustic Image Classification And Species Allocation) [START_REF] Petitgas | Sampling variance of species identification in fisheries-acoustic surveys based on automated procedures associating acoustic images and trawl hauls[END_REF] that directly infers species mixing proportions at the echogram level from global features computed as global statistics on the echogram. The method is based on hierarchical classification of echograms, where each class refers to a trawled site or a group of trawled sites. The fish abundance of a given image is estimated as the proportion of the class to which it belongs. Global methods are more appropriate than school based ones because trawled sites involve a fish mixture and the use of the echogram provides a direct relation between acoustic data and species proportions. Our proposed method is a generalization of the AICASA method. As AICASA, we use global feature estimated on fish schools extracted from echograms (typically, 3 nautical milelong echograms) using MOVIES software [START_REF] Weil | MOVIES-B: an acoustic detection description software: application to shoal species classification[END_REF]. But, unlike AICASA, we do not classify the echograms and attribute the same mixing proportions for all images belonging to the same class; we estimate the inference of species mixing proportions for each echogram as a weighted sum of the mixing proportions evaluated at trawled sites. The weights are computed according to a fast marching algorithm [START_REF] Sethian | Fast Marching Methods and Level Set Methods for Propagating Interfaces[END_REF] based on a distance between echogram statistics. Second, in AICASA, no spatial regularization is enforced, whereas some spatial regularity is expected, this is a built-in feature of the proposed approach. The proposed method is applied on synthetic and real data and evaluated in terms of mixing proportion estimation accuracy. The presented paper is organized as follows: in Section II, we introduce the problem and in Section III we describe the proposed method. In section IV, we give some results on synthetic and real data and in Section V we draw some conclusions.

II. PROBLEM STATEMENT

Backscattered-acoustic signal is digitized providing acoustic samples that form the echogram. The continuous recorded echogram is divided into bins. The binning results in a set of acoustic images that covers the survey area (typically, 3 nautical mile-long echograms). Figure 1 shows an example of echograms. be this area to which we associate a spatial grid (see figure 2). Each grid point s I corresponds to an acoustic image (points in black in figure 2). In some grid points (presented by circles in figure 2), denoted by { } K k k

I : 1 =
where K is the total number of trawl catches, we have information about fish abundance through catch hauls. Each catch indexed by k is mainly composed of a mixture of several species denoted by a vector k Π of J elements, where J is the number of fish species and each vector element

) ( j k Π
represents the proportion of the species j in the trawl catch k such that:

k j J j k ∀ = Π ∑ =1 , 1 ) ( (1)
∑ = Π = Π K k k s k s I w 1 ) ( (2)
III.

METHOD DESCRIPTION

The weights { } K k k w : 1 = estimation is based on the hypothesis that acoustic image features are associated to the fish species mixture in the image and on the fact that fish abundance has some consistency in space. We propose an interpolation method based on the model given by equation ( 2 and based also on the similarity of their acoustic content. Acoustic content of each image is given by a set of acoustic parameters computed from fish schools extracted from each image. School objects were identified and extracted using MOVIES software [START_REF] Sethian | Fast Marching Methods and Level Set Methods for Propagating Interfaces[END_REF]. From the extracted fish school we compute a set of feature distributions. As in [START_REF] Scalabrin | Narrow band acoustic identification of mono specific fish shoals[END_REF], we use parameters belonging to four categories: school positions, school morphology, school density and the occupation of space by schools. In this application we used a set of 20 parameters proposed in [START_REF] Scalabrin | Narrow band acoustic identification of mono specific fish shoals[END_REF]. Acoustic features are computed as the distribution of the different school features estimated on each image (mainly 100 schools). The distributions are estimated according to Parzen method [START_REF] Parzen | On the Estimation of a Probability Density Function and the Mode[END_REF] and are denoted by { }

Ω ∈ = s f s f I H , 20 : 1 ) (
where f is the feature index.

According to these acoustic features, we compute the similarity between two images as follows:

∑ = = 20 1 )) ( ), ( ( ) , ( f s f k f s k I H I H KL I I Dist ( 3
)
where KL is the Kullback-Leibler divergence [START_REF] Kullback | On information and sufficiency[END_REF].

We use this distance within a fast marching algorithm [START_REF] Sethian | Fast Marching Methods and Level Set Methods for Propagating Interfaces[END_REF] to

compute the weights { } K k k w : 1 =
. Fast marching algorithm is a numerical technique designed to track the evolution of interfaces. It aims at marching outward, computing distances from mesh points to the interface. The key step is to compute the distance to the next nearest mesh point. Then the front moves further outward with velocity F according to the Hamilton-Jacobi equation:

F T 1 = ∇ (4)
where 0 > F is the front moving speed and T is the travel time.

= = K i T i T k T k I T G I T G I w 1 )) ( ( )) ( ( ) ( (5)
where G is a radial kernel.

IV.

Results and discussion

We applied the proposed method on both synthetic and real data.

A. Synthetic data Synthetic data are generated from real fish school features as follows:

-We estimate fish features (e.g., fish density, target strength, depth, school shape as a set of ellipsoidal volumes); -We simulate several fish school using OASIS acoustic fish school simulator [START_REF] Diner | Correction on school geometry and density: approach based on acoustic image simulation[END_REF], which generates acoustic echograms from a given echosounder and vessel setting and the given school characteristics; -We fix a target fish proportion grid; -We associate to each grid point an echogram comprising schools randomly selected from the simulated dataset according to set class proportions.

The random selection induces a 'noisy' proportion map. During the test, we suppose known the proportions in some images selected randomly (identified images) and then, we apply the proposed interpolation method to estimate the abundance on the rest of images.

We show two grid examples that are simulated with real school dataset from Ifremer [START_REF] Scalabrin | Narrow band acoustic identification of mono specific fish shoals[END_REF]:

-a four-region case for which each region is composed of a single specie: 'Sardine', 'Anchovy', 'Horse Mackerel' and 'Chinchard' ; -a six-region case for which each region is a mixture of mainly two species. In figures 3 and 4, we show the initial data grids respectively a 128x128 and 250x250 squared images (each pixel refers to an acoustic image). For each identified image we show the fish abundance k Π . In figure 5, we show the simulated proportion map for the sardine specie for the six-region case. It might be noticed that the simulation procedure results in "noisy" proportion maps. In table I, we give the estimation error mean for each specie and for both cases. The error is computed as follows:

( )

∑ Ω ∈ Π - Π = T s e j s s j j I I ) ( ) ( , , τ (6)

V. Conclusion

We proposed a new for fish proportion estimation from acoustic features and some trawl catches. The interpolation is based on a weighting scheme computed from fast marching algorithm. The method allows taking into account both of the acoustic similarity between images and also the spatial distance between non identified echograms and the trawled areas.

In future work other global features extracted directly from the echogram and not only on the schools extracted from it, can be investigated such as the echogram backscatter distribution and texture.

Figure 1 .

 1 Figure 1. An echogram example (MOVIES software)

Figure 2 .

 2 Figure 2. Notations

  ) that computes the weight based on the spatial proximity of non-identified images { } s I to the identified ones { }

F 1 and

 1 T may also be interpreted as travel cost and intrinsic instance respectively. The time T at a grid point s I depends on the punctual velocity ) ( s I F and also on the spatial distance between the point s I position and the starting point position. We propose to use K fast marching processes, each initialized by an identified image k used as the weights for the fish abundance estimation.

∑

Figure 3 .

 3 Figure 3.4 monospecific case

Figure 4 .

 4 Figure 4 . 6bispecific case

Figure 5 .

 5 Figure 5. Sardine proportion map

Figure 6 .

 6 Figure 6 .Sardine estimation error map B. Real data Figures 7-10 show the interpolation results of the fish proportions ('Anchovy', 'Sardine', 'Horse Mackerel' and 'Sprat') for the French fisheries-acoustic survey of IFREMER performed on spring 2000 over the French shelf of Biscay performed with the RV "Thalassa". Trawl catch positions are presented in red circles.

Figure 7 .Figure 8 .Figure 9 .Figure 10 .

 78910 Figure 7. Anchovy proportion map

Table I .

 I Estimation error rate for each specie proportion and for each example.