Data-Driven Upscaling o f Orientation Kinematics i n Suspensions o f Rigid Fibres

Adrien Scheuer , Amine Ammar, Emmanuelle Abisset-Chavanne Elias Cueto,Francisco Chinesta, Roland Keunings and Suresh G. Advani

Abstract: 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 st ate-of-the-art closures proved to be inadequate or simply do not exist.

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

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 themicroscopic scale,theorientation of asingleparticleisidentifi ed by aunit vector p aligned with theparticleaxis.In thecaseof a Newtonian suspending fl uid,theevolution of therod orientation is governed by Jeffery’sequation[Jeffery(1922)]

whith ?v the unperturbed fluid velocity g radient. These kinematics, derived under the assumptions of a Stokes fl ow, 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, Chinesta et al. (2016)] 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 effi ciently 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, 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 fi bres, 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

wheretheintegration is performed on theunit sphere S.

Theseorientation tensorsexhibit particular properties.a and A arecompletely symmetric,that is

and

Moreover,dueto thenormalisation condition ofψ,thetraceof a equals1.

Thesecond-order orientation tensor hasan intuitivephysical interpretation.A high valueof the fi rst(resp.second and third)diagonal component of a indicates that the particles tend to orient along this x-(resp.y-and z-)direction.If all thediagonal components are,the 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 tothe tensor suggests two-dimensional random or planar biaxial orientation state.Finally,aunit diagonal component indicatesfull alignment in that direction.

Thetimeevolutionof thesecond-order orientation tensor isreadily 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. [Jack and Smith (2004, 2007)].

?Simple closures: linear (LIN) [Hand (1962)] (exact in the case of isotropic orientations), quadratic (QUAD) [Hinch and Leal (1976)] (exact for aligned fibres) and hybrid (HYBR) [Advani and Tucker (1987)] (combining the two previous ones);

?Composite closures: attempting to approximate directly the second-order tensor A : ?v[Hinch and Leal (1976)];

?Orthotropic closures: attempting to express A in the principal axis of a [Cintra and Tucker (1995)];

?Natural closures: natural (NAT) [Dupret and Verleye (1999)] and IBOF [Chung and Kwon (2002)] are fitted closures based on the most general expansions of A in terms of a and ?v;

?Neural-newtork-based closures: NNET [Jack, Schache and Smith (2010)] and NNORT [Qadir and Jack (2009)];

?Closures for the sixth-order orientation tensor: such as LIN6, QUAD6, HYBR6[Advani and Tucker (1990)], or even invariant-based fitted closures INV6and IBF6[Jack and Smith (2005, 2006)].

The macroscopic scale offers a simpleand crude description of themicrostructure.Simulations at that scale are thus much cheaper,explaining why this description is preferred in industrial applications.Thepdf issubstituted by someof itsmoments,sacrificing thelevel 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. [Peherstorfer and 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, 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: S ince 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:

?Offl ine stage:construction of the data-base.A large amount of microscopic simulations involving populations of N fi bres are run,exploring a wide range of initial orientation confi gurations.Moreover,the different scenarios may also include variationsin thesuspension parameters,such asthefl ow velocity gradient,thefi brevolume fraction(infl uencing the inter-particle interactions),the confi nement state(see below),the applied electric fi eld in the case of charged particles(see below);these parametersarecollectively referred to asα.For each population of fi bres,themacroscopic descriptors,(a,a˙,α),are computed from the microscopic ones,(pi,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 itemsin thedatabase[Gonzálea,Chinesta and Cueto(2018)].

Figure1:Data-driven approach to fi bre orientation kinematics

?Onlinestage:macroscopic data-driven simulation.At each timestep,we identify in thedatabasetheclosest itemsto thecurrent orientation state(and suspension parameters)a(t,α)and combinethem to obtain theinstantaneousevolution 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.

Our aim is to present a general framework for the upscaling of orientation kinematics in suspensionsof rigid fibres.Moreadvanced and sophisticated techniquescanbeused toperform effi ciently both offl ineand onlinestages.On theonehand,theoffl ineconstruction of thedatabasecan berecast in thegeneral context of Design of Experiments(DoE)methodology[Eriksson,Johansson,Kettaneh-Wold et al.(2000)],to help decidewhich scenarios of interest to consider.Latin hypercubesampling[McKay,Beckman and Conover(1979)]is a popular sampling technique to explore parameter spaces in experimental designs.On theother 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 thedatabaseentries.

3 Illustration of theframework

In thissection,we propose an illustration of the methodology that hasjust 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 fi rst discusstheclassical unconfi ned dilutecase,for which thereisalong history of macroscopic models,allowing us to assess the performance of our approach.We will then discuss the relevance of the method for confi ned suspensions,for which traditional macroscopic models fail.Finally,we will discuss the case of semi-concentrated suspensions of electrically-charged fi bres,using microscopic direct numerical simulations inspired by molecular dynamics.

3.1 Unconfi ned dilute suspensionsof rodsimmersed in a simpleshear fl ow

In the caseof dilutesuspensions of rods immersed in a Newtonian fl uid,thekinematics of each particle follow Jeffery’s equation,Eq.(1).In this section,we consider that the suspension undergoesasimpleshear fl ow,whosevelocity fi eld isgiven by v=?˙γz 0 0?T,with˙γ=1 s-1.In this fl ow,it is known that the fi bres simply align in the fl ow fi eld(x-direction).

3.1.1 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 confi gurations,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 detailsand the Matlab implementation.In thisstudy,weconstruct two databases,thefi rst with(nm,ns)=(10,15)(150 initial confi gurations),and the second,more comprehensive with(nm,ns)=(20,12)(240 initial configurations).The nmindividual mean values are uniformly distributed over thesphere and the nsvariances range from 0.05(fi bresnearly aligned)to 1.75(fi bresnearly uniformly distributed all over the sphere).For each population,we run the fl ow simulation during 30 seconds,and compute each 10-1second the macroscopic descriptors a and˙a from the individual piand˙pi(i=1,...,N)as

and

3.1.2 Data-driven simulation

Duringtheonlinestage,weidentify inthedatabasethe K items adbk(k=1,...,K)closest to the current orientation tensor a(t).

Figure2:Exampleof distributionsof initial orientationsused to build thedatabase

To do so, we use the euclidean distance applied to the vectorized forms of the orientation

tensors, composed of its independent components:

Then,wecomputeaweighted averageof the K corresponding kinematics a˙dkbto derivethe instantaneousevolution a˙ that can beapplied to havetheorientation tensor at thenext time step,that is

Thereconstruction weights areobtained by solving theminimization problem

such that

Remark:The choice of the definition of distance is definitely a delicate q u estion. 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).

Figure 3: 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)

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).

3.1.3 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 a11of the second-order orientation tensor over ncrandom initial distributions composed of N fibres. The average L2 relative error is defined as

with

where amacro(t)is computed either using the data-driven approach or a closure-based model.As before,we consider the QUAD,HYBR and IBOFclosure approximations and amicro(t)is obtained form the expensive discrete microscopic simulations.The error estimator hereis based on the fi rst diagonal component of the orientation tensor only,sincein thecaseof ashear fl ow,thesuspending fi brestend to align in thefl ow fi eld,and thus a11is an approximateindicator of the alignment state asit approachesunity.

The results of this comparative study(nc=500)are shown in Fig.4.We observe that thedata-driven approach shows improved performancecompared 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(nm,ns)=(20,12))improvestheaccuracy of thedata-driven method.

3.2 Confi ned dilutesuspensionsof rodsimmersed in a simpleshear fl ow

Wenow movetoconfi ned suspensionsof rods,thatissuspensionsfl owingingapsnarrower than the fi brelength.Thegap wallsnow prevent theparticlefrom rotating freely and some trajectory,passing outside the fl ow domain are thus forbidden.We have shown in previous work[Perez,Scheuer,Abisset-Chavanne et al.(2016);Scheuer,Abisset-Chavanne,Chinesta et al.(2016,2018)]that in that case,the fi bre kinematics can be written as Jeffery’s equation augmented with an additional term that prevents the fi bre from leaving the fl ow domain,

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

Database constructionThe 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 simulationThe 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, 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.

Figure 5: 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)

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 assessmentWe use the same method as before to assess the performances of the method (here nc= 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).

3.3 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-, y- and z-directions, uniaxial elongation in the x- and y-directions and rotation flow around the x-, y- and 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.

3.4 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,Barinsiskiet al.(2015)].A microscopic model governing the evolution˙p of asinglerod isobtained from amicromechanical derivation using adumbbell representation of the particle and reads˙p=˙pJ+˙pE,where˙pEdepends on the external electric fi eld E and the charge q of the rod dipole.The proposed macroscopic model is however tainted with non-reliable closureapproximationsthat make it impractical to use.Wecould of courseapply thesamerationaleagain,building on top of thismodified Jeffery equation.However,we want to consider here a semi-concentrated suspension,that is we want to account for the effects of fi bre-fi bre interactions as well.Moreover,we want to emphasizethat thedata-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 simulationsinspired by molecular dynamics(MD).

This fi ne scale simulation technique is based on the following assumptions: (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 VLJand two other potentials, VEand VB, 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 confi guration of the particles.As before,we consider here a simpleshear fl ow,whosevelocity field isgiven by v=?˙γz 0 0?T,with˙γ=1 s-1.The electric fi eld 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 fi eld,the fi bres tend to align in the fl ow fi eld,as illustrated in Fig.7b,that shows the fi nal orientation state of the fi bres when E=0 NC-1.Conversely,when the electric fi eld is strong,the fi bres cannot align in the fl ow and the fi nal orientation is along an inclined axis(the inclination depends on the intensity of theelectric fi eld),asillustrated in Fig.7b,that showsthefi nal orientation stateof thefi bres when E=50 NC-1.

Database constructionIn 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 variationof theshear rate˙γor theinitial orientationstate(asintwopreviousillustrations)is astraightforward extension.Therefore,thedatabasesarebuiltby followingtheevolutionof populations of N=100,200,...800 particles subjected to an electric fi eld E of intensity ranging from 0 to 60 NC-1.Theshear rateisfi xed,˙γ=1 s-1and theisotropic orientation state is always chosen as initial configuration.In each case,we run the MD simulation during 10 seconds,and compute each 10-1second the macroscopic descriptors a and˙a from the individual piand˙pi(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.

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

Data-driven simulationThe 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: T he number of particles in the system N (that influences the amount of fibre-fibre interactions) and the intensity of th e e x ternal 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 (Nt) and the current value of the electric field intensity (Et). Among the databases at our disposal, we then identify the ones that best match the value of the current parameters (for example if (Nt, Et) = (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 w ay to handle for e xample 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 assessmentAgain, 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 nc= 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.

4 Conclusion and p erspectives

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

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-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, Chinesta et al. (2016)]). 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.

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