Learning the Macroscopic Flow Model of Short Fiber Suspensions from Fine-Scale Simulated Data

Fiber–fiber interaction plays an important role in the evolution of fiber orientation in semi-concentrated suspensions. Flow induced orientation in short-fiber reinforced composites determines the anisotropic properties of manufactured parts and consequently their performances. In the case of dilute suspensions, the orientation evolution can be accurately described by using the Jeffery model; however, as soon as the fiber concentration increases, fiber–fiber interactions cannot be ignored anymore and the final orientation state strongly depends on the modeling of those interactions. First modeling frameworks described these interactions from a diffusion mechanism; however, it was necessary to consider richer descriptions (anisotropic diffusion, etc.) to address experimental observations. Even if different proposals were considered, none of them seem general and accurate enough. In this paper we do not address a new proposal of a fiber interaction model, but a data-driven methodology able to enrich existing models from data, that in our case comes from a direct numerical simulation of well resolved microscopic physics.


Introduction
Short fiber reinforced composite is being used in a wide range of sectors in our society. The composite is manufactured from processes such as, among many others, injection molding, compression molding and extrusion. The mechanical performances of the manufactured composite, such as stiffness and strength, is determined by the final state of fiber orientation. The final fiber state results from the flow induced orientation evolution that in the case of semi-concentrated or concentrated flow regimes involves strong fiber-fiber interactions whose modeling is crucial for determining the final part properties and performance.
Predicting flow induced anisotropy in reinforced polymer forming processes is of major interest. However, most of the existing models, as the ones described below, have important limitations associated with the numerous hypotheses introduced during their analytical derivation. When proceeding out of the domain of validity of those models, a pragmatic approach consists of assuming the existing models are still valid and calibrating them from experimental tests.
When these calibrated nominal models are unable to reproduce experimental findings, some ad hoc terms are usually added, being then calibrated accordingly (e.g., diffusion term reflecting fiber interactions).
Thus, data was traditionally employed for calibrating models derived from physical considerations. However, many times those models failed to address experimental findings even when they were finely calibrated. It is the case of models describing intense fiber interactions in semi-concentrated and concentrated fiber suspensions. This issue motivated numerous works referred to later, and remains even today not totally solved.
Recent advances in artificial intelligence, and more particularly in machine learning, are opening new routes. Thus, data can be used not only for calibrating pre-assumed models, i.e., assimilating data into existing models to identify the most adequate values of their parameters. Data can be also used for inferring the model itself, or eventually the correction of a nominal pre-existing model to enhance its predictability capacity.
The present paper aims at addressing this important and timely topic of data-driven model enrichment, here addressed from a methodological viewpoint. Future works will focus on (i) the use of data-driven models for improving computational efficiency in the simulation of complex flows as the ones encountered in processing; (ii) the derivation of more accurate rheological models for describing complex fluids; and (iii) the use of these more accurate rheological models in processing simulations for circumventing, or at least alleviating, the limitation of current models and procedures and their associated predictions.

Revisiting Standard Flow-Induced Orientation Models
The orientation description and its flow induced evolution can be addressed at three different scales [1]:

•
At the microscopic scale the orientation state is described from the orientation of each fiber. Thus, in the case of a population consisting of N fibers of length 2L with an almost infinite aspect ratio (length to diameter), all of them are located into a representative volume element -RVEω, the orientation state in ω at time t reduces to the knowledge of the N unit vectors p i (t), centered at the i-fiber centre of gravity G i and aligned along its axis.
If we assume that population immersed in a fluid flow characterized by the velocity gradient ∇v, assumed unperturbed by the fibers' presence and orientation, when neglecting the fiber-fiber interactions, the fiber rotary velocity is expressed by the Jeffery equation [2], which readṡ • The mesoscopic description makes use of the probability distribution function Ψ(x, t, p), giving the fraction of fibers that at position x and time t are oriented along p.
Again, when neglecting fiber interaction, the evolution of that distribution function is given by the Fokker-Planck equationΨ whereΨ represents the material derivative and ∇ p (•) the gradient operator involving the conformational coordinates p.

•
The main difficulty of the previous descriptions is their complexity, the microscopic one due to the fact of needing to track all the fibers immersed in the flow and the mesoscopic due to the fact that the distribution function is defined in a high-dimensional space involving the physical space x, the time t and the orientation p, the last defined on the unit sphere. Macroscopic descriptions circumvent both difficulties by operating with the probability distribution function moments [3].
Thus, the second moment of the probability distribution function, also known as second order orientation tensor a(x, t), reads a = p ⊗ pΨ dp, where the integral applies on the unit sphere surface and "⊗" denotes the tensor product.
As described in [1], the second order moment time evolution, again ignoring fiber interaction, readsȧ where A is the fourth order moment and ":" represents the tensor product twice contracted.
In what follows the orientation rate in Equation (4) will be referred to asȧ Je f f for indicating that it is associated to the Jeffery model that ignores fiber-fiber interactions.
The macroscopic description is retained in most of works and commercial software because of its relative simplicity. However, for employing it, one must consider closure relations [4] for expressing the fourth order orientation tensor as a function of the second one, whose validity remains nowadays an open issue. In the numerical tests performed in the present paper the IBOF closure [5] (invariant-based orthotropic fitted) is considered. Moreover, by restricting the orientation description to the second moment, only some probability distribution functions can be accurately described, others cannot be expressed from the solely use of the second order orientation tensor. On the other hand, size effects were neglected when deriving the models at different scales, limiting the model validity when addressing confined flows.
In the discussion above, fiber-fiber interaction was ignored. Pioneering works assumed that interactions can be considered as randomizing events that tends to push the systems towards an isotropic orientation state, and then, they were modeled from a diffusion term into the Fokker-Planck equation by Folgar and Tucker [6]. Thus, the probability distribution function time evolution is now expressed fromΨ with D the diffusion coefficient, assumed scaling the flow intensity from the second invariant of the velocity gradient. By using this expression it results in an orientation evolution equation composed of two terms, the usual Jeffery-based orientation model complemented with a correction that can be easily derived as illustrated in [1]. Thus the time evolution of the orientation can be written in its most general form aṡ where the first term corresponds to the Jeffery contribution given by Equation (4) and the second contribution, the deviation with respect to the Jeffery model, comes from the interactions, modeled as a diffusion term within the Folgar and Tucker (F&T) model. For addressing anomalous experimental orientation findings at odds with F&T predictions, in particular a delay in the orientation rates, other models of fiber-fiber interactions have been proposed [7][8][9][10][11][12]. In these studies, the observed delay was attributed to the intense fiber-fiber interactions and modeled either with a modified diffusion term, a microscopic description of interactions, or by introducing a sliding mechanism between fibers and fluid, with an anisotropic rotary diffusion.
At present neither, closure relations nor interaction models that are general and accurate enough, exist. A recent comparison of different approaches can be found in [13]. The accuracy of the fiber orientation is tested usually by predicting the effect on mechanical properties due to the error in predicting of fiber orientation. This issue is beyond the scope of this work but references [3,4] have done that comparison. In [14] authors noticed that for example, the incorrect prediction of fiber orientation as the suspension flows into a rib could seriously compromise the part performance.
Despite the fact that numerous works continue to introduce improved models, the present paper proposes an alternative methodology for improving existing models, based on the data-driven enriching of the model.
In our former works [15] we addressed the fully data-driven modeling of suspensions, here we restrict to enriching the Jeffery model with a data-driven description of interactions, a methodology also considered in some of our former works [16,17].
As our main aim is not the derivation of a new interaction model for fiber suspensions, in what follows we will consider a molecular dynamics (MD) simulation, able to represent in absence of interactions Jeffery motion, and in presence of interactions describing them with a a standard Lennard-Jones potential. Here we are only interested in assimilating that data extracted from the data-driven interaction model that when added to the Jeffery contribution enables an accurate prediction of the orientation evolution in flow conditions that are different from the ones used to extract the model. This fine-scale simulated data-generator does not pretend to infer an improved model for simulating processes, but only proving that a macroscopic orientation model can be learned from microscopic data in order to accurately describe that orientation state. By having access to more accurate data associated with the orientation evolution of a population of rods, from adequate measurements or generated from more accurate microscopic models (direct numerical simulations solving the fluid kinematics and induced fibers motion) more accurate macroscopic models could be derived.
The methodology described in this paper remains totally general, the only difference being the quality of the data that serves for learning the model, and consequently it can be applied to any flow regime, from dilute (without major interest being the Jeffery model valid) to the highly concentrated regime (where standard model fail).
The next section addresses the data-generator based on a MD fine-scale simulation. Then in Section 3 we introduce the machine learning techniques that will be applied in Section 4 to extract the interaction term. Finally in Section 5 we proceed to verify the proposal by conducting some numerical experiments.

Data Generator Based on Molecular Dynamics (MD) Simulations
The fine scale simulation is performed by representing the fiber as M connected beads where forces apply, the last derived from different potentials as well as the hydrodynamic drag.
We consider three different potentials: The Lennard-Jones potential V LJ to describe bead-bead interaction (both beads belonging to different fibers) and V E and V B that describe the elongation and bending interaction between neighbor beads of the same fiber.
As indicated at the end of the previous section, the MD simulation here addressed allows recovering Jeffery model predictions in absence of interaction. When interactions occur nothing was demonstrated on the ability of MD for representing fibers motion, but here as previously indicated, our main aim is proving that a macroscopic model can be extracted from individual microscopic data. The quality of the learned model only depends on the quality of the data that served to extract it. The MD strategy here considered must be viewed as a simply synthetic data generator.
Thus, without loss of generality, in what follows we consider the following potential for addressing the just mentioned mechanisms where N is the number of fibers in the population and V LJ i,j;k,l represents the Lennard-Jones potential between bead i of fiber j and bead k of fiber l. The characteristic function χ j,l = 1 − δ jl (with δ the Kroenecker delta) ensures that the Lennard-Jones potential only acts between beads belonging to different fibers.
The Lennard-Jones potential reads with d ik = ||x k − x i || the distance between beads B i of fiber j and bead B k of fiber l, located at positions x i and x k respectively, and and σ are the two usual parameters involved in the Lennard-Jones potential. The elongation potential is given by where d eq is the equilibrium distance between two successive beads, distance at which the potential reaches its minimum (and consequently the associated force vanishes). In the previous expression K E reflects the potential intensity. The bending potential between three successive beads where θ i−1,i,i+1 is the angle defined by vectors joining beads i − 1, i and i, i + 1 with θ eq the angle at equilibrium, that vanishes when considering straight fibers. In our study we consider parameters ensuring an almost rigid behavior of the fiber. The hydrodynamic drag was enforced by adding at each bead the external force whereẋ i denotes the velocity of bead B i and ξ the friction coefficient. The total force applied on each bead, inserted into the Newton equation, allows the acceleration calculation, and from it the velocity and position update, as performed in standard MD simulations.
The problem is expressed in a dimensionless form prior to proceed with its solution. It is important to emphasize that in absence of Lennard-Jones interactions and for almost rigid fibers, the previous modeling framework recovers the usual Jeffery kinematics.

Methods
Two dimensionality reduction methods, the principal component analysis (PCA) and the Code2Vect will be employed for constructing the interaction data-driven model, the former for reducing the dimensionality and the last for creating the nonlinear regression. For the sake of completeness both are summarized in the next sections.

Principal Component Analysis (PCA)
Let us consider a vector y ∈ R D containing some data. If some correlations exists in that data, one could find a linear transformation W to transform y into x, with x ∈ R d , with d < D. This mapping reads y = Wx.
The transformation matrix W, D × d, satisfies the orthogonality condition W T W = I d , where I d represents the d × d-identity matrix (WW T is not necessarily I D ).
Assume that there exist M different snapshots y 1 , . . . , y M , which we store in the columns of a D × M matrix Y. The associated d × M reduced matrix X contains the associated vectors x i , i = 1, . . . , M. We assume without loss of generality centered variables.
PCA proceeds by guaranteeing maximal preserved variance and enforcing the latent variables in x being uncorrelated. In other words, the covariance will be diagonal. PCA extracts the d uncorrelated latent variables from Pre-and post-multiplying by W T and W, respectively, and making use of the fact that W T W = I d , gives us The covariance matrix C yy can then be factorized by applying the singular value decomposition (SVD), i.e., C yy = VΛV T , with V containing the orthonormal eigenvectors and Λ the diagonal matrix containing the eigenvalues, sorted in descending order. Substituting Equation (16) into Equation (15), we arrive at with C xx diagonal as soon as the d columns of W are taken collinear with d columns of V.
In the data-driven modeling approach proposed later, PCA will be used for expressing a time history in a time reduced basis. Thus, instead of describing a time function from its value at different time steps, it will be described by the weight of few functions defined in the whole time interval involved in the time reduced basis.

Code2Vect
In [18] the authors proposed a supervised classification and nonlinear regression procedure that was called Code2Vect because of its ability to map data from a representation (non-metric) space into a vector space equipped of a metric.
We assume that points in the origin space consists of M arrays composed on D entries, noted by y i . Theirs images in the vector space are noted by x i ∈ R d , with in general d D. The mapping is described by the d × D matrix W ensuring where both, the components of W and the images x i ∈ R d , i = 1, · · · , M, must be calculated. Each point x i keeps the label (value of the quantity of interest (QoI)) associated with is origin point y i , denoted by O i . We would like placing points x i , such that the Euclidian distance with each other point x j scales with their outputs difference, i.e., where the coordinates of one of the points x i can be arbitrarily chosen. Linear mappings are limited and do not allows proceeding in nonlinear settings. Thus, a better choice consists of the nonlinear mapping W(y), expressible from the general polynomial form where W k are d × D matrices and P (y) a polynomial basis. Equation (19) ensures clustering due to the fact that data with similar output will be placed closer than when the outputs differ significantly. As well as the mapping has been determined, Equation (18) (or its nonlinear counterpart) becomes an efficient nonlinear regression as described in [18].
In the data-driven modeling approach proposed later, Code2Vect will be used for defining nonlinear regressions eventually involving heterogeneous data (e.g., the components of the orientation and velocity gradient tensors).

Model Training
In what follows the orientation based on the Jeffery model that ignores interaction will be denoted by the superscript • Je f f , whereas the orientation evaluated from MD simulations will be denoted by • MD .
The deviation • dev is expected to describe the fiber interactions because the proposed MD model in absence of interactions (with a vanishing Lennard-Jones potential) was able to reproduce very accurately fiber motion as described by Jeffery's equation (Equation (4)) . The associated rate equation can be written asȧ It is expected that in absence of confinement, size effect, and for a given fiber concentration,ȧ dev will depend on the flow velocity gradient and the orientation state a itself.
In what follows we consider two scenarios: • The first is related to the orientation process from a given initial orientation state for different flows characterized by their velocity gradient kept constant during all the flow process; • The second is more general and looks for expressing the deviation rate as a function of both the applied flow (defined by its changing velocity gradient) and the existing orientation state at a given time instant.

Orientation Evolution in Homogeneous Flows
In this case, the deviation rate at any time is assumed depending on the initial orientation state a 0 and the applied flow ∇v. Our main purpose is determining accurately a(t; a 0 , ∇v). For the sake of simplicity and without loss of generality, we assume the same initial orientation, almost isotropic.
The procedure for learning the model consists of the following steps: • Consider an almost initial isotropic fiber distribution p k (t = 0), k = 1, . . . , N, defining the initial orientation tensor a 0 ≈ I/3.

•
Generate randomly or using a design of experiments (DoE) different flows characterized by the velocity gradients ∇v j , j = 1, . . . , M.
• The second invariant of tensors ∇v j ,γ j , allows defining the maximum simulation time for each flow t j max such thatγ j t j max = Γ, with Γ total strain large enough for reaching an almost steady orientation state.
• Perform the MD integration for all the considered flows from the initial fiber distribution p k (t = 0), leading to p j k (t), from which the second order orientation tensor a MD j (t) can be computed at time t from that enables calculating its time derivativeȧ DM j (t).
• Perform the integration of the Jeffery model from a 0 for the different flows, to obtain a Je f f j (t) and its associate orientation rateȧ The interval Γ is divided in D intervals, and the deviation rate at the associated deformation (the deformation is defined from γ =γt) extracted, i.e.,ȧ dev j (γ l ), l = 1, . . . , D.

•
Thus, each component (m, n) ofȧ dev j (γ l ),ȧ dev mn,j (γ l ), can be expressed as a point in R D , y dev mn,j , whose l−component resultsȧ dev mn,j (γ l ). • Now, the PCA is applied and the first d modes retained. For facilitating visualization we will consider in fact the first three modes, i.e., d = 3, defining the matrix W. Thus, the image of y dev mn,j is given by z dev mn,j = Wy dev mn,j , with z dev mn,j ∈ R d . • Now, Code2Vec applies on vectors G j (with G j the vector form of the velocity gradient ∇v j ), that are transformed into g j from the nonlinear mapping W mn (G j ), i.e., where the mapping is constructed by enforcing With orientation model W mn (G) constructed, as soon as a new orientation trajectory must be evaluated for a given velocity gradient G it suffices to proceed as follows: • Map G into the vector space according to g = W(G)G; • Identify the set S containing the closest neighbors g i of g; • Compute the approximation weights ω i , i = 1, . . . , S, depending on the distance between g i and g, for example by using radial basis; • Interpolate the deviation reduced orientation trajectories according to • Push up to reconstruct the time evolution according to It is important to note that the inverse mapping R d → R D is not defined, and by this reason we perform the same interpolation that in the reduced space (previous item)

Instantaneous Orientation Behavior
The present case is even simpler. We consider different orientation states, that can be computed by taking some snapshots of the fibers orientation at different time during a complex flow simulation. Then, for each of them, we applied different flows to evaluate the instantaneous orientation rate of change.
The procedure reads: • Define at t = 0, Q initial orientation states p q k , k = 1, . . . , N & q = 1, . . . , Q, using for example a MD simulation in a complex flow. • Compute the associated orientation tensors a q 0 • Define J flows defined by their velocity gradient tensor ∇v j , j = 1, . . . , J. Now, online, for any a and ∇v, it results y that by applying the mapping leads to z mn = W mn y. As in the previous section, the set of closest neighbors can be defined, and an interpolation performed to obtainȧ dev , that is used in the time-marching integration procedure according tȯ that adds to the standard Jeffery contribution the one extracted from the data related to the deviations.

Results
In the considered suspension the fiber volume fraction (φ) is kept at 0.01 (semi-dilute regime, 0.0025 < φ < 0.05 for a fiber aspect ratio of 20). These conditions apply to both procedures previously defined.

Orientation Evolution in Homogeneous Flows
Total 35 cases with different ∇v j are prepared for this case study (some of them related to pure shear and elongation flows). Out of 35 cases, its aimed to predict theȧ dev of the case number 10. The initial orientation state is almost isotropic, whereas the velocity gradient is The prediction for the 10th flow is made according to the procedure described in Section 3.3.1 and is shown in Figure 1. The mapped point representing the 10th flow has as closest neighbors flows 4, 7, 9 and 11. The velocity gradient ∇v j for j = 4, 7, 9&11 are: In what follows, and for the sake of representation clarity, only the component • 11 of the orientation tensor will be discussed.
The correspondingȧ dev 11 (t) for the five flows j = 4, 7, 9, 10, 11 are presented in Figure 2a, whereas the prediction and reference solution (original data) for that flow, j = 10, are presented in Figure 2b. The proposed strategy perform quite reasonably by capturing the main behavior. Obviously the finer is the sampling, the more accurate is the prediction.

Instantaneous Orientation Behavior
In the present case we consider two numerical tests reported below.
• In this first test, we perform a MD simulation for extracting different orientation states when considering an almost isotropic initial orientation state a the flow characterized by The extracted Q = 4 intermediate orientation states are defined in Table 1.  For these four orientation states a total of 38 samples at different ∇v are prepared.
As discussed in Section 3.3.2 Code2Vect applies on the M data ∇v j and a q 0 whose results are presented in Figure 3, where original data and predictions ofȧ dev 11 up toγt = 2 is depicted (left). Prediction at the evaluation point (star) along with 37 samples (filled circle) is depicted on the right image of that figure, where color scales with the value ofȧ dev 11 atγt = 0.3. The predicted data shows a good agreement with the raw data validating the proposed methodology.

•
The last numerical test addresses a more complex scenario where the velocity gradient changes within the simulation-time window.
The prediction forȧ dev is made when fibers (almost isotropically distributed at the initial time) are placed under two velocity gradient fields, ∇v 1 (shear) and ∇v 2 (elongation), consecutively, the former applying during the first 3.3 s and the last the subsequent 10 s.
The velocity gradient ∇v for these flow fields is Data-driven enriched Jeffery solution and the one provided by MD simulation are compared over the entire time period in Figure 4a,b, whereas Figure 4c depicts the deviation prediction (following the procedure described in Section 3.3.2) versus the reference one provided directly by the data. As it can be noticed, the numerical values of the prediction when compared with the reference data, both reported in Table 2, exhibit an excellent accuracy.

Conclusions and Prospects
This work shows the opportunity of hybrid modeling methodologies, combining a model based on physics complemented with a data-driven enrichment to account for the gap between the model predictions and the collected data. In this work the synthetic data was used from a high-fidelity fine-scale simulation based on molecular dynamics, a simulation which took fiber-fiber interactions into account.
The data-driven model was constructed by using advanced machine learning techniques able to operate at the low-data limit, in particular the PCA and the so-called Code2Vect.
The proposed methodology was verified with numerical experiments that show the efficiency of our hybrid approach and the techniques employed for constructing the data-driven model.
Currently injection molding simulations use a constitutive model in which they have to provide an interaction parameter which is the same everywhere and is determined from experiments conducted under simple shear. It has been shown that the flow field is very different around corners, inserts, ribs and undercuts in complex injection molded parts. It is expected that the constitutive interaction parameter locally will be different but there is no model or methodology to determine this. Our approach should allow for more accurate fiber orientation prediction in such situations. The extension of the present work to efficient processing simulations constitute a work in progress.