Data-Driven Upscaling of Orientation Kinematics in Suspensions of Rigid Fibres

Describing the orientation state of the particles is often critical in fibre suspension applications. Macroscopic descriptors, the so-called second-order orientation tensor (or moment) leading the way, are often preferred due to their low computational cost. Closure problems however arise when evolution equations for the moments are derived from the orientation distribution functions and the impact of the chosen closure is often unpredictable. In this work, our aim is to provide macroscopic simulations of orientation that are cheap, accurate and closure-free. To this end, we propose an innovative data-based approach to the upscaling of orientation kinematics in the context of fibre suspensions. Since the physics at the microscopic scale can be modelled reasonably enough, the idea is to conduct accurate offline direct numerical simulations at that scale and to 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 an imprecise macroscopic model. This methodology is presented in the well-known case of dilute fibre suspensions (where it can be compared against closure-based macroscopic models) and in the case of suspensions of confined or electrically-charged fibres, for which state-of-the-art closures proved to be inadequate or simply do not exist.


Data-Driven Upscaling of Orientation Kinematics in Suspensions of Rigid Fibres 1 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 [Jeffery (1922)] p J = ∇v · p − (∇v : (p ⊗ p))p, 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 [Chinesta (2013)], bending phenomena [Abisset-Chavanne, Férec, Ausias et al. (2015)], particle inertia [Scheuer, Grégoire, Abisset-Chavanne et al. (2018)], electrical forces [Perez, Abisset-Chavanne, Barinsiski et al. (2015)], wall effects [Perez, Scheuer, Abisset-Chavanne et al. (2016); Scheuer, Abisset-Chavanne, ] is readily achievable using a dumbbell representation of a rod [Bird, Curtiss, Armstrong et al. (1987); Binetruy, Chinesta and Keunings (2015)]. 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 [Lozinski and Chauvière (2003); Chauvière and Lozinski (2004)]. Notable progress were made recently [Chinesta, Ammar, Leygue et al. (2011)] with the introduction of the Proper Generalized Decomposition [Ammar, Mokdad, Chinesta et al. (2006], 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. [Advani and Tucker (1987)], read respectively and where the integration is performed on the unit sphere S. These orientation tensors exhibit particular properties. a and A are completely symmetric, that is and 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) 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. Smith (2004, 2007)].
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) [Darema (2004); Michopoulos, Farhat and Houstis (2004)], entailing the ability to incorporate additional data into an executing application, improving modelling methods and prediction capabilities. Kirchdoerfer et al. [Kirchdoerfer and Ortiz (2016)] 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. Willcox (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 [Perez, Scheuer, Abisset-Chavanne et al. (2016); Scheuer, Abisset-Chavanne, ]. 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, ȧ , α), are computed from the microscopic ones, (p i , ṗ 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 ȧ = f (a, α) from interpolations of the items in the database [Gonzálea, ].  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 [Eriksson, Johansson, Kettaneh-Wold et al. (2000)], to help decide which scenarios of interest to consider. Latin hypercube sampling [McKay, Beckman and Conover (1979)] 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.

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 [Chen (2015)] 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

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 .
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 The reconstruction weights are obtained by solving the minimization problem 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 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 [Perez, Scheuer, Abisset-Chavanne et al. (2016); Scheuer, Abisset-Chavanne, ] 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, p =ṗ 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 [Perez, Scheuer, Abisset-Chavanne et al. (2016); Scheuer, Abisset-Chavanne, ], 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. 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 [Perez, Abisset-Chavanne, Barinsiski et al. (2015)]. 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. [Perez, Scheuer, Abisset-Chavanne et al. (2018)]. 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-  Figure 8: 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 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 [Brunn (1977)] for second-order fluid in the limit of low Weissenberg and the recent multi-scale modelling based on that model [Borzacchiello, Abisset-Chavanne, ). 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.