Direct simulation of fiber suspensions under smalland large amplitude oscillatory shear flow

Rabih Mezher1*

College of Engineering and Technology, American University of the Middle East, Kuwait

*Corresponding author: Rabih Mezher, College of Engineering and Technology, American University of the Middle East, Kuwait. E-mail: Rabih.Mezher@aum.edu.kw

Citation: Rabih Mezher1 (2020) Direct simulation of fiber suspensions under small and large amplitude oscillatory shear flow. Int J Mater Res Sci Tech 1(1): 25-34.

https://dx.doi.org/10.47890/IJMRST/2020/RMezher/14205434

Received Date: 24 July, 2020; Accepted Date: 30 July, 2020; Published Date: 3 August , 2020

Abstract

Experimentally, oscillatory shear flow is used to investigate polymeric viscoelastic properties. However, the knowledge of the microstructure’s orientation evolution (kinematics) with time is essential in order to predict those properties, and is difficult experimentally. Therefore it is often required to conduct numerical simulations to get insight on the kinematics evolution of fiber suspensions. The main aim of the current work is to use direct numerical simulations to investigate the time evolution of the orientation and kinematics of fiber suspensions under small and large amplitude oscillatory shear flow.

Keywords: fiber suspensions; direct numerical simulation; oscillatory shear flow; fiber-fiber interactions;

Introduction

Today, fiber suspensions are important materials in a wide range of applications. In medical applications, for example, nanocomposites from natural fibers are materials that can be used in cardiovascular implants, in articular cartilage and in tissue engineering [1]. In the composite industry, fiber suspensions take a major part due to their interesting mechanical, electrical, thermo-physical properties and their cost efficient processing [2]. Indeed, the majority of thermoplastic and polymeric materials have poor mechanical properties and are often reinforced with fibers. Depending on the complexity of the material, local variations can appear in the fiber orientation which can affect the mechanical, thermal, and insulative properties. It is therefore recommended that one be able to predict the flow behavior of the composite fluid in connection to the fiber orientation in order to optimize mold design and the performance of the material [3, 4]. Moreover, since the bulk rheological properties of fiber suspensions depend on the filler and the matrix, the knowledge of these properties is also important to deter- mine the internal structure. In general, the rheological behavior of the filled systems is quite complex. It is difficult to predict these behaviors experimentally, therefore simulating fiber suspensions immersed in a flow has been the main goal of several previous works and publications [5–10].

In general, the theory of fiber suspensions consists of two main components: an equation describing the evolution of the fiber orientation and an equation related to the suspensions rheology relating the stress and the rate of strain. There are three concentration regimes of fiber suspensions: diluted, semi-diluted and concentrated [11]. When the concentration is low enough (diluted regime), the kinematic of each fiber is given from the fluid kinematics unperturbed by the presence of the fibers and moreover, it does not depend on the existence of all the others. In this case one can describe the suspensions by tracking a population of fibers that moves with the suspending fluid and orients depending on the velocity gradient according to the Jeffery’s equation [12] that relates the orientation evolution with the flow velocity field. Thus each fiber evolves independently of the others in terms of kinematics and orientation evolution. However when the concentration increases, the difficulty to describe such systems rises. In the semi-diluted or semi-concentrated regimes, the fibers rarely touch but are subjected to long-range hydrodynamic interactions from the perturbation of fluid motion caused by nearby fibers [13, 14]. And in this case, Folgar and Tucker [15] added a term to Jeffery’s equation representing an isotropic diffusivity proportional to the effective deformation rate to account for fiber–fiber interactions influencing the degree of fiber alignment. When the concentration is high enough (concentrated regimes) fiber–fiber interactions can no longer be neglected due to the large number of interactions occurring between the fibers and appropriate models describing these interactions must be considered. A number of studies focused on fibers immersed in a Newtonian fluid where fiber-fiber interactions were neglected [16, 2]. Thus, those works cannot describe the effect of intense interactions between the fibers in a concentrated regime and they are only valid for diluted or semi-diluted suspensions, which is not the case of many industrial processes such as compression of Sheet Molding Compounds or Glass Mat Transfer or injection of Bulk Molding Compounds. Moreover, in the concentrated regime models, operating on the macroscopic (allowing fast simulations) do not allow to address the fine physics involved in such intense interactions. Few works focused on proposing microscopic descriptions in order to better describe complex fiber–fiber inter- actions [17, 18].

One way of studying the kinematics of fiber suspensions and their macroscopic properties is to conduct a direct numerical simulation (or DNS for abbreviation). A DNS is based on the computation, in a representative volume, of the motion of hundreds of fibers and their interactions. It is a step by step process which derives kinematics as well as macroscopic properties, while taking into account the forces applied on each fiber [19]. DNS allows a precise description of the fine physics involved at the microscopic scale.

Polymeric viscoelastic properties are commonly investigated through the use of oscillatory shear experiments in order to have insight on the microstructure [20]. There are many studies in the literature done on the rheological behavior of filled polymers, which exhibit remarkable non-linear viscoelasticity even at very low deformation amplitude. In general, the addition of filler increases the storage and loss moduli of the neat resin [21–24], especially when the fiber aspect ratio (length/diameter) is large [20]. Few works used DNS in order to better describe the complex kinematics of fiber suspensions and fiber-fiber interactions, under oscillatory shear flows [17, 18, 25, 26].

The objective of the current work is to use DNS to investigate the time evolution of the orientation and kinematics of fiber suspensions under small and large amplitude oscillatory shear flow (SAOS and LAOS). In the following sections, the micro-mechanical model based on DNS is thoroughly presented. Dilute and concentrated suspensions are considered. Unlike previous works [27–30] where fibers had the same aspect ratio, some simulations will consider fibers having different lengths. In fact, the assumption of fibers with different lengths is close to real industrial composites [24]. Then the time evolution of the kinematics of the fibers immersed in a Newtonian fluid under small and large amplitude oscillatory shear flow is studied.

Micro-mechanical model

Main assumptions of the Mirco-mechanical model

The main assumptions and hypothesis of the DNS model are:

  • Dilute and concentrated fiber suspensions are considered;
  • The suspending fluid is Newtonian, incompressible and the flow is laminar;
  • The velocity field and its gradient are assumed to be homogeneous in the considered volume, and unperturbed by the fluid’s presence;
  • The mass of the fibers is negligible, thus the inertia of the fibers are not considered in this study;
  • Fibers are considered to be prolate spheroids in shape and have the same diameter;
  • Interactions between fibers are considered through contact and lubrication forces;
  • Initially and before the simulation is started, the fibers are uniformly dis- tributed in the considered volume and they are not in any state of inter- action. Their orientation state is initially almost isotropic;
  • In this model, each fiber (denoted by α) have a length l(α) and a diameter d. (Figure 1). Thus the aspect ratio of the fibers can be defined by r(α) such as:
r α = l α d                  (1) MathType@MTEF@5@5@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGYbWdamaaCaaaleqabaWdbmaabmaapaqaa8qacqaHXoqyaiaa wIcacaGLPaaaaaGccqGH9aqpdaWcaaWdaeaapeGaamiBa8aadaahaa Wcbeqaa8qadaqadaWdaeaapeGaeqySdegacaGLOaGaayzkaaaaaaGc paqaa8qacaWGKbaaaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaae iiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqG GaGaaeiiaiaabccacaqGGaGaaeikaiaabgdacaqGPaaaaa@4E51@

Figure 1: Fiber modeled by a prolate spheroid

When the suspensions are considered dilute, then the volume concentration φ is related with the aspect ratio as:

ϕ 1 r 2                (2) MathType@MTEF@5@5@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacqaHvpGzcqGHKjYOdaWcaaWdaeaapeGaaGymaaWdaeaapeGaamOC a8aadaahaaWcbeqaa8qacaaIYaaaaaaakiaabccacaqGGaGaaeiiai aabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGa aeiiaiaabccacaqGGaGaaeiiaiaabIcacaqGYaGaaeykaaaa@483C@

Where r is the average aspect ratio. When the suspensions are considered concentrated, then the volume concentration φ is related with the aspect ratio as:

ϕ 1 r                (3) MathType@MTEF@5@5@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacqaHvpGzcqGHLjYSdaWcaaWdaeaapeGaaGymaaWdaeaapeGaamOC aaaacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaae iiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqG OaGaae4maiaabMcaaaa@473C@

The more ris important (i.e. long fibers), the more the system is considered concentrated for a fixed value of φ. In this model the fibers are supposed such as 1>r.

Considerations for fiber position andorientation

The position of the center of gravity Gfor a fiber denoted by α is at r(α) (Figure 2) with:

r α = x α e x + y α e y + z α e z                (4) MathType@MTEF@5@5@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGYbWdamaaCaaaleqabaWdbmaabmaapaqaa8qacqaHXoqyaiaa wIcacaGLPaaaaaGccqGH9aqpcaWG4bWdamaaCaaaleqabaWdbmaabm aapaqaa8qacqaHXoqyaiaawIcacaGLPaaaaaGccaWGLbWdamaaBaaa leaapeGaamiEaaWdaeqaaOWdbiabgUcaRiaadMhapaWaaWbaaSqabe aapeWaaeWaa8aabaWdbiabeg7aHbGaayjkaiaawMcaaaaakiaadwga paWaaSbaaSqaa8qacaWG5baapaqabaGcpeGaey4kaSIaamOEa8aada ahaaWcbeqaa8qadaqadaWdaeaapeGaeqySdegacaGLOaGaayzkaaaa aOGaamyza8aadaWgaaWcbaWdbiaadQhaa8aabeaak8qacaqGGaGaae iiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqG GaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGOaGaaeinaiaabM caaaa@5DF2@

ex, ey and ez represent respectively the three unit vectors related to the three space coordinates (in the three dimensional axis system ). x(α), y(α) and z(α) are the coordinates of the center of gravity of fiberα. The direction or orientation of a fiber is defined by the unit vector p(α) such as p α = p 1 α x+ p 2 α y+ p 3 α z MathType@MTEF@5@5@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGWbWdamaaCaaaleqabaWdbmaabmaapaqaa8qacqaHXoqyaiaa wIcacaGLPaaaaaGccqGH9aqpcaWGWbWdamaaDaaaleaapeGaaGymaa WdaeaapeWaaeWaa8aabaWdbiabeg7aHbGaayjkaiaawMcaaaaakiaa dIhacqGHRaWkcaWGWbWdamaaDaaaleaapeGaaGOmaaWdaeaapeWaae Waa8aabaWdbiabeg7aHbGaayjkaiaawMcaaaaakiaadMhacqGHRaWk caWGWbWdamaaDaaaleaapeGaaG4maaWdaeaapeWaaeWaa8aabaWdbi abeg7aHbGaayjkaiaawMcaaaaakiaadQhaaaa@50B1@ p1(α), p2(α) and p3(α)

Are the coordinates of the orientation vector. For a population of fibers, a statistical orientation distribution function Ψ(p) describes the state of the fibers orientation in the fluid. In the literature, the orientation state of a population of fibers is given by the second order orientation tensor a2 [31].

a ik = p p i p k ψ p dp              (5) MathType@MTEF@5@5@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaqGHbWdamaaBaaaleaapeGaamyAaiaadUgaa8aabeaak8qacqGH 9aqpdaqfqaqabSWdaeaapeGaamiCaaqab0WdaeaapeGaey4kIipaaO GaaeydGiaadchapaWaaSbaaSqaa8qacaWGPbaapaqabaGcpeGaamiC a8aadaWgaaWcbaWdbiaadUgaa8aabeaak8qacqaHipqEdaqadaWdae aapeGaamiCaaGaayjkaiaawMcaaiaadsgacaWGWbGaaeiiaiaabcca caqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiai aabccacaqGGaGaaeiiaiaabccacaqGOaGaaeynaiaabMcaaaa@5440@

Where p is the unit vector along the fiber’s axis characterizing its orientation, while pi and pk are two components of p; aik are the components of the second order orientation tensor.

Figure 2: Schematic representation of the fiber in the spatial system of the laboratory

Index i and k denote the direction of space (i=1, 2 or 3 and k=1, 2 or 3). Orientation tensors give a statistical average for the orientation of a fiber population [31]. The orientation tensor, preserves a unit trace through time, and is fully symmetric (i.e.aik=aki). In the 3 dimensional (3D) case it is written explicitly in the matrix form:

a 2 = a 11 a 12 a 13 a 21 a 22 a 23 a 31 a 32 a 33                       (6)      MathType@MTEF@5@5@+= feaahqart1ev3aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWHHbWdamaaBaaaleaapeGaaGOmaaWdaeqaaOWdbiabg2da9maa dmaapaqaauaabaqadmaaaeaapeGaaeyya8aadaWgaaWcbaWdbiaaig dacaaIXaaapaqabaaakeaapeGaaeyya8aadaWgaaWcbaWdbiaaigda caaIYaaapaqabaaakeaapeGaaeyya8aadaWgaaWcbaWdbiaaigdaca aIZaaapaqabaaakeaapeGaaeyya8aadaWgaaWcbaWdbiaaikdacaaI XaaapaqabaaakeaapeGaaeyya8aadaWgaaWcbaWdbiaaikdacaaIYa aapaqabaaakeaapeGaaeyya8aadaWgaaWcbaWdbiaaikdacaaIZaaa paqabaaakeaapeGaaeyya8aadaWgaaWcbaWdbiaaiodacaaIXaaapa qabaaakeaapeGaaeyya8aadaWgaaWcbaWdbiaaiodacaaIYaaapaqa baaakeaapeGaaeyya8aadaWgaaWcbaWdbiaaiodacaaIZaaapaqaba aaaaGcpeGaay5waiaaw2faaiaabccacaqGGaGaaeiiaiaabccacaqG GaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabc cacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeii aiaabccacaqGOaGaaeOnaiaabMcacaqGGaGaaeiiaiaabccacaqGGa Gaaeiiaaaa@67EC@

Since a discrete population of a large number of fibers is present is this work, the integral in expression (5) can be replaced by a discrete sum. Hence, for a given quantity of n fibers this tensor is calculated as:

a ik = 1 n α=1 n p i α p k α                    (7) MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaqGHbWdamaaBaaaleaapeGaamyAaiaadUgaa8aabeaak8qacqGH 9aqpdaWcaaWdaeaapeGaaGymaaWdaeaapeGaamOBaaaadaGfWbqabS WdaeaapeGaeqySdeMaeyypa0JaaGymaaWdaeaapeGaamOBaaqdpaqa a8qacqGHris5aaGccaqGnaIaamiCa8aadaqhaaWcbaWdbiaadMgaa8 aabaWdbmaabmaapaqaa8qacqaHXoqyaiaawIcacaGLPaaaaaGccaWG WbWdamaaDaaaleaapeGaam4AaaWdaeaapeWaaeWaa8aabaWdbiabeg 7aHbGaayjkaiaawMcaaaaakiaabccacaqGGaGaaeiiaiaabccacaqG GaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabc cacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGOaGaae4n aiaabMcaaaa@5D69@

A fixed axis system (0,x,y,z) is defined relative to the reference of the laboratory and a local axis system (G,x’,y’,z’)is defined relative to the reference of each fiber (Figure 2).

Representative elementary volume

DNS is done here in a representative elementary volume (or REV for abbreviation) containing the fibers and the fluid. The volume considered is a cubic reference cell (Figure 3). The main aim is to reproduce the bulk behavior of the composite; but since it is not possible to simulate a very large number of fibers, this difficulty can be overcomed by replicating the cubic cell throughout space at each time step [32, 33, 34]. Consequently if a fiber moves in the original cell, its image in each of the neighboring cells moves in exaclty the same way.

Figure 3: Representative elementary volume

It follows that if a fiber leaves the reference cell, one of its images will enter through the opposite face.

The fiber in the original cell can interact with any of the fibers in the surrounding cells (Figure 3 where fibers interacts with the green image fibers) . The positions of the fibers are updated at each time step only in the original cell that is then replicated.

Fiber motion: translation

The velocity field of the suspending fluid, is imposed by considering an oscillatory shear flow where the shear γ and the shear rate γ are given as a function of time t by [35]:

γ= γ 0 sin wt , γ ˙ = γ 0 wcos wt                   (8) MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacqaHZoWzcqGH9aqpcqaHZoWzpaWaaSbaaSqaa8qacaaIWaaapaqa baGcpeGaci4CaiaacMgacaGGUbWaaeWaa8aabaWdbiaadEhacaWG0b aacaGLOaGaayzkaaGaaiilaiqbeo7aN9aagaGaa8qacqGH9aqpcqaH ZoWzpaWaaSbaaSqaa8qacaaIWaaapaqabaGcpeGaam4Daiaadogaca WGVbGaam4Camaabmaapaqaa8qacaWG3bGaamiDaaGaayjkaiaawMca aiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGa GaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabcca caqGGaGaaeiiaiaabIcacaqG4aGaaeykaaaa@5D59@

for an amplitude γ0 and w=2πf (f is the frequency of oscillation). When the amplitude of the oscillation (γ0) is small enough, the stress response of the fluid remains sinusoidal, and the fluid remains in its linear regime. At these amplitudes the stress response is linear with the amplitude [24]. These tests are called small-amplitude oscillatory shear (SAOS). If the amplitude is increased, the stress response of a fluid may no longer be sinusoidal. For these large-amplitude oscillatory shear (LAOS) tests, a nonlinear response may be obtained for the stress [24].

The velocity field is:

V T x = V 1 , V 2 , V 3 = γ ˙ y,0,0                  (9)   MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGwbWdamaaCaaaleqabaWdbiaadsfaaaGcdaqadaWdaeaapeGa amiEaaGaayjkaiaawMcaaiabg2da9maabmaapaqaa8qacaWGwbWdam aaBaaaleaapeGaaGymaaWdaeqaaOWdbiaacYcacaWGwbWdamaaBaaa leaapeGaaGOmaaWdaeqaaOWdbiaacYcacaWGwbWdamaaBaaaleaape GaaG4maaWdaeqaaaGcpeGaayjkaiaawMcaaiabg2da9maabmaapaqa a8qacuaHZoWzpaGbaiaapeGaamyEaiaacYcacaaIWaGaaiilaiaaic daaiaawIcacaGLPaaacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaa bccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaae iiaiaabccacaqGGaGaaeiiaiaabIcacaqG5aGaaeykaiaabccacaqG Gaaaaa@5B6C@

(throughout the entire text, “T” denotes the transpose of any vector or matrix) and is defined in each point of space.

Axis x is associated with the first direction of the flow velocity field (V1 in equation (9)), axis y is associated with the velocity gradient and axis z is associated with the vorticity (Figure 2). The velocity gradient is given by:

V= V m x n = 0 γ ˙ 0 0 0 0 0 0 0                   (10) MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacqGHhis0caWGwbGaeyypa0ZaaSaaa8aabaWdbiabgkGi2kaadAfa paWaaSbaaSqaa8qacaWGTbaapaqabaaakeaapeGaeyOaIyRaamiEa8 aadaWgaaWcbaWdbiaad6gaa8aabeaaaaGcpeGaeyypa0ZaamWaa8aa baqbaeaabmWaaaqaa8qacaaIWaaapaqaa8qacuaHZoWzpaGbaiaaae aapeGaaGimaaWdaeaapeGaaGimaaWdaeaapeGaaGimaaWdaeaapeGa aGimaaWdaeaapeGaaGimaaWdaeaapeGaaGimaaWdaeaapeGaaGimaa aaaiaawUfacaGLDbaacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaa bccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaae iiaiaabccacaqGGaGaaeiiaiaabccacaqGOaGaaeymaiaabcdacaqG Paaaaa@5AF3@

m=1, 2 and 3 defines the index for the three components of V(x), and n=1, 2 and 3 defines the index for the three space coordinates. When the velocity field is uniform in the considered volume, then the flow velocity field at position r(α) can be written as:

V r α =V r α                  (11)  MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGwbWaaeWaa8aabaWdbiaadkhapaWaaWbaaSqabeaapeWaaeWa a8aabaWdbiabeg7aHbGaayjkaiaawMcaaaaaaOGaayjkaiaawMcaai abg2da9iabgEGirlaadAfacqGHflY1caWGYbWdamaaCaaaleqabaWd bmaabmaapaqaa8qacqaHXoqyaiaawIcacaGLPaaaaaGccaqGGaGaae iiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqG GaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabI cacaqGXaGaaeymaiaabMcacaqGGaaaaa@55A4@

The rate of strain tensor is written as:

D= 1 2 V+ (V) T                  (12)     MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGebGaeyypa0ZaaSaaa8aabaWdbiaaigdaa8aabaWdbiaaikda aaWaaeWaa8aabaWdbiabgEGirlaadAfacqGHRaWkcaGGOaGaey4bIe TaamOvaiaacMcapaWaaWbaaSqabeaapeGaamivaaaaaOGaayjkaiaa wMcaaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccaca qGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaa bccacaqGGaGaaeikaiaabgdacaqGYaGaaeykaiaabccacaqGGaGaae iiaiaabccaaaa@539A@

and the vorticity tensor is given by:

W= 1 2 V (V) T                  (13) MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGxbGaeyypa0ZaaSaaa8aabaWdbiaaigdaa8aabaWdbiaaikda aaWaaeWaa8aabaWdbiabgEGirlaadAfacqGHsislcaGGOaGaey4bIe TaamOvaiaacMcapaWaaWbaaSqabeaapeGaamivaaaaaOGaayjkaiaa wMcaaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccaca qGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaa bccacaqGGaGaaeikaiaabgdacaqGZaGaaeykaaaa@512D@

The position of G relative to the fluid is given by a vector that will be denoted by q. The relative velocity of G with respect to the fluid is thus given by:

q ˙ α = r ˙ α V r α = r ˙ α γ ˙ y α x                (14)   MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qaceWGXbWdayaacaWaaWbaaSqabeaapeWaaeWaa8aabaWdbiabeg7a HbGaayjkaiaawMcaaaaakiabg2da9iqadkhapaGbaiaadaahaaWcbe qaa8qadaqadaWdaeaapeGaeqySdegacaGLOaGaayzkaaaaaOGaeyOe I0IaamOvamaabmaapaqaa8qacaWGYbWdamaaCaaaleqabaWdbmaabm aapaqaa8qacqaHXoqyaiaawIcacaGLPaaaaaaakiaawIcacaGLPaaa cqGH9aqpceWGYbWdayaacaWaaWbaaSqabeaapeWaaeWaa8aabaWdbi abeg7aHbGaayjkaiaawMcaaaaakiabgkHiTiqbeo7aN9aagaGaa8qa caWG5bWdamaaCaaaleqabaWdbmaabmaapaqaa8qacqaHXoqyaiaawI cacaGLPaaaaaGccaWG4bGaaeiiaiaabccacaqGGaGaaeiiaiaabcca caqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiai aabccacaqGGaGaaeiiaiaabIcacaqGXaGaaeinaiaabMcacaqGGaGa aeiiaaaa@6485@

q(α) is associated with the resultant force acting on the fiber (for the detailed expressions of r(α) and q(α) the reader can refer to [19], [36] and [37]). Assuming that the forces of inertia in translation and rotation are negligible with respect to the hydrodynamic forces, the equation of the fiber’s motion is reduced to:

F α + ζ α q ˙ α =0                (15)     MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGgbWdamaaCaaaleqabaWdbmaabmaapaqaa8qacqaHXoqyaiaa wIcacaGLPaaaaaGccqGHRaWkcqaH2oGEpaWaaWbaaSqabeaapeWaae Waa8aabaWdbiabeg7aHbGaayjkaiaawMcaaaaakiabgwSixlqadgha paGbaiaadaahaaWcbeqaa8qadaqadaWdaeaapeGaeqySdegacaGLOa GaayzkaaaaaOGaeyypa0JaaGimaiaabccacaqGGaGaaeiiaiaabcca caqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiai aabccacaqGGaGaaeiiaiaabccacaqGOaGaaeymaiaabwdacaqGPaGa aeiiaiaabccacaqGGaGaaeiiaaaa@58DC@

Where F(α) is the resultant force acting of the fiber and ς(α) is the resistance tensor relating the applied force F(α) and the relative velocity q(α) in the shear flow. When they are no short-range forces acting on the fiber, it would move in the flow like a force-free single independent fiber (i.e. according to the Jeffery’s equation; for further details, the reader can refer to [19]). Ωis the rotating velocity of the fluid in the flow field (Vr(α)) given by:

Ω= 1 2 ×V = 1 2 ε ijk x j V k                (16)  MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWHPoGaeyypa0ZaaSaaa8aabaWdbiaaigdaa8aabaWdbiaaikda aaWaaeWaa8aabaWdbiabgEGirlabgEna0kaadAfaaiaawIcacaGLPa aacqGH9aqpdaWcaaWdaeaapeGaaGymaaWdaeaapeGaaGOmaaaacqaH 1oqzpaWaaSbaaSqaa8qacaWGPbGaamOAaiaadUgaa8aabeaak8qada WcaaWdaeaapeGaeyOaIylapaqaa8qacqGHciITcaWG4bWdamaaBaaa leaapeGaamOAaaWdaeqaaaaak8qacaWGwbWdamaaBaaaleaapeGaam 4AaaWdaeqaaOWdbiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeii aiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGa GaaeiiaiaabIcacaqGXaGaaeOnaiaabMcacaqGGaaaaa@5C93@

Where “×” is the cross product and εijk is the Levi-Civita permutation tensor.

Interaction between two fibers

The configuration of interaction is defined by the distance or gap between the fibers, when interactions happen. For two fibers α and β, the normal vector orthogonal to their plane through the orientations p(α) and p(β) is written as:

n α,β = p α × p β p α × p β              (17)  MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGUbWdamaaCaaaleqabaWdbmaabmaapaqaa8qacqaHXoqycaGG SaGaeqOSdigacaGLOaGaayzkaaaaaOGaeyypa0ZaaSaaa8aabaWdbi aadchapaWaaWbaaSqabeaapeWaaeWaa8aabaWdbiabeg7aHbGaayjk aiaawMcaaaaakiabgEna0kaadchapaWaaWbaaSqabeaapeWaaeWaa8 aabaWdbiabek7aIbGaayjkaiaawMcaaaaaaOWdaeaapeWaaqWaa8aa baWdbiaaygW7aiaawEa7caGLiWoacaWGWbWdamaaCaaaleqabaWdbm aabmaapaqaa8qacqaHXoqyaiaawIcacaGLPaaaaaGccaaMb8Uaey41 aqRaaGzaVlaadchapaWaaWbaaSqabeaapeWaaeWaa8aabaWdbiabek 7aIbGaayjkaiaawMcaaaaakmaaemaapaqaa8qacaaMb8oacaGLhWUa ayjcSdaaaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabc cacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGOaGaaeym aiaabEdacaqGPaGaaeiiaaaa@6D1B@

Figure 4: Contact or lubrication between two fibers α and β

The distance from the mass center of fiber β to the interaction point with fiber is l(β,α) (Figure 4). In the case of high aspect ratio prolate spheroids, the gap between two fibers Θ(α,β) (Figure 4) is [38]:

Θ α,β = r α r β n α,β d 2 14 l (α,β) 2 l (α) 2 d 2 14 l (β,α) 2 l (β) 2             (18)   MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaqGyoWdamaaCaaaleqabaWdbmaabmaapaqaa8qacqaHXoqycaGG SaGaeqOSdigacaGLOaGaayzkaaaaaOGaeyypa0ZaaeWaa8aabaWdbi aadkhapaWaaWbaaSqabeaapeWaaeWaa8aabaWdbiabeg7aHbGaayjk aiaawMcaaaaakiabgkHiTiaadkhapaWaaWbaaSqabeaapeWaaeWaa8 aabaWdbiabek7aIbGaayjkaiaawMcaaaaaaOGaayjkaiaawMcaaiab gwSixlaad6gapaWaaWbaaSqabeaapeWaaeWaa8aabaWdbiabeg7aHj aacYcacqaHYoGyaiaawIcacaGLPaaaaaGccqGHsisldaWcaaWdaeaa peGaamizaaWdaeaapeGaaGOmaaaadaGcaaWdaeaapeGaaGymaiabgk HiTiaaisdadaWcaaWdaeaapeGaamiBa8aadaahaaWcbeqaa8qacaGG OaGaeqySdeMaaiilaiabek7aIjaacMcapaWaaWbaaWqabeaapeGaaG Omaaaaaaaak8aabaWdbiaadYgapaWaaWbaaSqabeaapeGaaiikaiab eg7aHjaacMcapaWaaWbaaWqabeaapeGaaGOmaaaaaaaaaaWcbeaaki abgkHiTmaalaaapaqaa8qacaWGKbaapaqaa8qacaaIYaaaamaakaaa paqaa8qacaaIXaGaeyOeI0IaaGinamaalaaapaqaa8qacaWGSbWdam aaCaaaleqabaWdbiaacIcacqaHYoGycaGGSaGaeqySdeMaaiyka8aa daahaaadbeqaa8qacaaIYaaaaaaaaOWdaeaapeGaamiBa8aadaahaa Wcbeqaa8qacaGGOaGaeqOSdiMaaiyka8aadaahaaadbeqaa8qacaaI YaaaaaaaaaaaleqaaOGaaeiiaiaabccacaqGGaGaaeiiaiaabccaca qGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeikaiaa bgdacaqG4aGaaeykaiaabccacaqGGaaaaa@841F@

Interaction forces

The resulting interaction force on fiber α is:

F α = μα F c α,μ n α,μ + βα F lb α,β n α,β             (19)  MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGgbWdamaaCaaaleqabaWdbmaabmaapaqaa8qacqaHXoqyaiaa wIcacaGLPaaaaaGccqGH9aqpdaGfqbqabSWdaeaaa8qabeqdpaqaa8 qacqGHris5aaGccaqGnaYaaSbaaSqaaiabeY7aTjabgcMi5kabeg7a HbqabaGccaWGgbWdamaaDaaaleaapeGaam4yaaWdaeaapeWaaeWaa8 aabaWdbiabeg7aHjaacYcacqaH8oqBaiaawIcacaGLPaaaaaGccaWG UbWdamaaCaaaleqabaWdbmaabmaapaqaa8qacqaHXoqycaGGSaGaeq iVd0gacaGLOaGaayzkaaaaaOGaey4kaSYaaybuaeqal8aabaaapeqa b0WdaeaapeGaeyyeIuoaaOGaaeydGmaaBaaaleaacqaHYoGycqGHGj sUcqaHXoqyaeqaaOGaaeydGiaadAeapaWaa0baaSqaa8qacaWGSbGa amOyaaWdaeaapeWaaeWaa8aabaWdbiabeg7aHjaacYcacqaHYoGyai aawIcacaGLPaaaaaGccaWGUbWdamaaCaaaleqabaWdbmaabmaapaqa a8qacqaHXoqycaGGSaGaeqOSdigacaGLOaGaayzkaaaaaOGaaeiiai aabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGa aeiiaiaabccacaqGGaGaaeikaiaabgdacaqG5aGaaeykaiaabccaaa a@7752@

Equation (19) shows that F(α) is divided into two sums: the first one groups all fibers denoted by μ in contact with fiber α and here Fc(α,μ) is the magnitude of the contact force; the second one groups all fibers denoted by β in lubrication with fiber and here Flb(α,β) is the magnitude of the lubrication force.

First interaction: lubrication forces

When two fibers come close together, a lubrication force occurs due to the fluid being trapped in between the fibers. The lubrication force Flb(α,β) was calculated within a linear approach given by Yamane et al [38]:

F lb α,β =K Θ α,β Θ α,β n α,β = F lb α,β n α,β                   (20)     MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGgbWdamaaDaaaleaapeGaamiBaiaadkgaa8aabaWdbmaabmaa paqaa8qacqaHXoqycaGGSaGaeqOSdigacaGLOaGaayzkaaaaaOGaey ypa0JaeyOeI0Iaam4samaabmaapaqaa8qacaqGyoWdamaaCaaaleqa baWdbmaabmaapaqaa8qacqaHXoqycaGGSaGaeqOSdigacaGLOaGaay zkaaaaaaGccaGLOaGaayzkaaGaaeiMd8aadaahaaWcbeqaa8qadaqa daWdaeaapeGaeqySdeMaaiilaiabek7aIbGaayjkaiaawMcaaaaaki aad6gapaWaaWbaaSqabeaapeWaaeWaa8aabaWdbiabeg7aHjaacYca cqaHYoGyaiaawIcacaGLPaaaaaGccqGH9aqpcaWGgbWdamaaDaaale aapeGaamiBaiaadkgaa8aabaWdbmaabmaapaqaa8qacqaHXoqycaGG SaGaeqOSdigacaGLOaGaayzkaaaaaOGaamOBa8aadaahaaWcbeqaa8 qadaqadaWdaeaapeGaeqySdeMaaiilaiabek7aIbGaayjkaiaawMca aaaakiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccaca qGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaa bccacaqGGaGaaeiiaiaabIcacaqGYaGaaeimaiaabMcacaqGGaGaae iiaiaabccacaqGGaaaaa@797F@

with:

K= 3π η 0 d 2 Θ α,β p α × p β          (21)   MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGlbGaeyypa0ZaaSaaa8aabaWdbiaaiodacqaHapaCcqaH3oaA paWaaSbaaSqaa8qacaaIWaaapaqabaGcpeGaamiza8aadaahaaWcbe qaa8qacaaIYaaaaaGcpaqaa8qacaqGyoWdamaaCaaaleqabaWdbmaa bmaapaqaa8qacqaHXoqycaGGSaGaeqOSdigacaGLOaGaayzkaaaaaO WaaqWaa8aabaWdbiaaygW7aiaawEa7caGLiWoacaWGWbWdamaaCaaa leqabaWdbmaabmaapaqaa8qacqaHXoqyaiaawIcacaGLPaaaaaGcca aMb8Uaey41aqRaaGzaVlaadchapaWaaWbaaSqabeaapeWaaeWaa8aa baWdbiabek7aIbGaayjkaiaawMcaaaaakmaaemaapaqaa8qacaaMb8 oacaGLhWUaayjcSdaaaiaabccacaqGGaGaaeiiaiaabccacaqGGaGa aeiiaiaabccacaqGGaGaaeiiaiaabIcacaqGYaGaaeymaiaabMcaca qGGaGaaeiiaaaa@683A@

where ηο is the viscosity of the suspending fluid, and Θ(α,β) is the relative velocity between two nearest interaction points, projected on the normal n(α,β) (for the detailed calculation of the lubrication forces the reader can refer to [39])

Second interaction: contact forces

Contact forces are assumed to occur if the gap between two close fibers is equal to zero and if Fc(α,β)≠0. In previous works, authors focused on introducing simple conditions to impose the contact, because a complex description of contact forces is diffcult to generate [29, 40]. The condition considered here uses equation (18) and takes the following form [19]:

d dt r α r β n α,β = Θ α,β =0                 (22)   MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qadaWcaaWdaeaapeGaamizaaWdaeaapeGaamizaiaadshaaaWaamWa a8aabaWdbmaabmaapaqaa8qacaWGYbWdamaaCaaaleqabaWdbmaabm aapaqaa8qacqaHXoqyaiaawIcacaGLPaaaaaGccqGHsislcaWGYbWd amaaCaaaleqabaWdbmaabmaapaqaa8qacqaHYoGyaiaawIcacaGLPa aaaaaakiaawIcacaGLPaaacqGHflY1caWGUbWdamaaCaaaleqabaWd bmaabmaapaqaa8qacqaHXoqycaGGSaGaeqOSdigacaGLOaGaayzkaa aaaaGccaGLBbGaayzxaaGaeyypa0JaaeiMd8aadaahaaWcbeqaa8qa daqadaWdaeaapeGaeqySdeMaaiilaiabek7aIbGaayjkaiaawMcaaa aakiabg2da9iaaicdacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaa bccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaae iiaiaabccacaqGGaGaaeiiaiaabIcacaqGYaGaaeOmaiaabMcacaqG GaGaaeiiaaaa@68D3@

with Θ(α,β)≈0. This condition is only valid if the fiber that is in contact has a null velocity. It physically means that two fibers in contact cannot overlap or pass through each other. They maintain a steady null distance between them at the contact point, when the gap is almost equal to zero (for the detailed calculation of the contact forces the reader can refer to [39]).

Fiber motion: rotation

The evolution of the fiber’s orientation in time is given by the equation:

p ˙ α = p α × ω α Ω        (23)   MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qaceWGWbWdayaacaWaaWbaaSqabeaapeWaaeWaa8aabaWdbiabeg7a HbGaayjkaiaawMcaaaaakiabg2da9iabgkHiTiaadchapaWaaWbaaS qabeaapeWaaeWaa8aabaWdbiabeg7aHbGaayjkaiaawMcaaaaakiab gEna0oaabmaapaqaa8qacqaHjpWDpaWaaWbaaSqabeaapeWaaeWaa8 aabaWdbiabeg7aHbGaayjkaiaawMcaaaaakiabgkHiTiaahM6aaiaa wIcacaGLPaaacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccaca qGGaGaaeikaiaabkdacaqGZaGaaeykaiaabccacaqGGaaaaa@54FB@

where Ω is the rotating velocity of the fluid in the flow field V(x) and ω(α) is the relative rotating velocity of the fiber with respect to that of the fluid. p(α) is written as:

p ˙ α = p ˙ J α + p ˙ I α       (24)    MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qaceWGWbWdayaacaWaaWbaaSqabeaapeWaaeWaa8aabaWdbiabeg7a HbGaayjkaiaawMcaaaaakiabg2da9iqadchapaGbaiaadaqhaaWcba WdbiaadQeaa8aabaWdbmaabmaapaqaa8qacqaHXoqyaiaawIcacaGL PaaaaaGccqGHRaWkceWGWbWdayaacaWaa0baaSqaa8qacaWGjbaapa qaa8qadaqadaWdaeaapeGaeqySdegacaGLOaGaayzkaaaaaOGaaeii aiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeikaiaabkdacaqG0a GaaeykaiaabccacaqGGaGaaeiiaaaa@5025@

Where pJ(α) is the Jeffery orientation evolution and pI(α) is the orientation evolution due to the presence of interactions (for the detailed calculation of p(α) and expression of (24) the reader can refer to [19], [36] and [37]). This relative rotation produces a torque written as ξ(α)(α) Another torque is caused by the deformation of the fluid due to the presence of the velocity field. This one is given by H α :D   MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaWefv3ySLgznf gDOfdaryqr1ngBPrginfgDObYtUvgaiuqaqaaaaaaaaaWdbiab=Tqi i9aadaahaaWcbeqaa8qadaqadaWdaeaapeGaeqySdegacaGLOaGaay zkaaaaaOGaaiOoaiaadseacaqGGaGaaeiiaaaa@4701@ . The two new resistance tensors ξ(α) and H α MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaWefv3ySLgznf gDOfdaryqr1ngBPrginfgDObYtUvgaiuqaqaaaaaaaaaWdbiab=Tqi i9aadaahaaWcbeqaa8qadaqadaWdaeaapeGaeqySdegacaGLOaGaay zkaaaaaaaa@442A@ , called the rotating resistance tensors represent the friction due to the rotation of the fibers and the deformation of the fluid respectively (for the detailed expressions of these tensors the reader can refer to [19] and [36]). Finally the resultant torque of the lubrication and contact forces, occurs when two fibers are in interaction (Figure 4). Denoting these two fibers by α and β, the resultant torque on one fiber is:

T α = l α,β p α × F α,β          (25)      MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGubWdamaaCaaaleqabaWdbmaabmaapaqaa8qacqaHXoqyaiaa wIcacaGLPaaaaaGccqGH9aqpcaWGSbWdamaaCaaaleqabaWdbmaabm aapaqaa8qacqaHXoqycaGGSaGaeqOSdigacaGLOaGaayzkaaaaaOGa amiCa8aadaahaaWcbeqaa8qadaqadaWdaeaapeGaeqySdegacaGLOa GaayzkaaaaaOGaey41aqRaamOra8aadaahaaWcbeqaa8qadaqadaWd aeaapeGaeqySdeMaaiilaiabek7aIbGaayjkaiaawMcaaaaakiaabc cacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeii aiaabIcacaqGYaGaaeynaiaabMcacaqGGaGaaeiiaiaabccacaqGGa Gaaeiiaaaa@5B7E@

where l(α,β) is the distance from the mass center of fiber to the interaction point with fiber β (Figure 4). Where F(α,β) is the resultant force (equation (19)) applied to the couple of fibers in interaction αand β. Neglecting the inertia of fibers, the balance of moments gives:

T α + ξ α ω α + H α :D=0             (26)    MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGubWdamaaCaaaleqabaWdbmaabmaapaqaa8qacqaHXoqyaiaa wIcacaGLPaaaaaGccqGHRaWkcqaH+oaEpaWaaWbaaSqabeaapeWaae Waa8aabaWdbiabeg7aHbGaayjkaiaawMcaaaaakiabgwSixlabeM8a 39aadaahaaWcbeqaa8qadaqadaWdaeaapeGaeqySdegacaGLOaGaay zkaaaaaOGaey4kaSYefv3ySLgznfgDOfdaryqr1ngBPrginfgDObYt UvgaiuqacqWFlecspaWaaWbaaSqabeaapeWaaeWaa8aabaWdbiabeg 7aHbGaayjkaiaawMcaaaaakiaacQdacaWGebGaeyypa0JaaGimaiaa bccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaae iiaiaabccacaqGGaGaaeiiaiaabccacaqGOaGaaeOmaiaabAdacaqG PaGaaeiiaiaabccacaqGGaaaaa@67BE@

From equation (26) one can find the total relative rotating velocity of the fiber with respect to the fluid ω(α):

ω α = ω J α + ω I α           (27)   MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacqaHjpWDpaWaaWbaaSqabeaapeWaaeWaa8aabaWdbiabeg7aHbGa ayjkaiaawMcaaaaakiabg2da9iabeM8a39aadaqhaaWcbaWdbiaadQ eaa8aabaWdbmaabmaapaqaa8qacqaHXoqyaiaawIcacaGLPaaaaaGc cqGHRaWkcqaHjpWDpaWaa0baaSqaa8qacaWGjbaapaqaa8qadaqada WdaeaapeGaeqySdegacaGLOaGaayzkaaaaaOGaaeiiaiaabccacaqG GaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabI cacaqGYaGaae4naiaabMcacaqGGaGaaeiiaaaa@547E@

where ω1(α) is the Jeffery relative rotating velocity and ω1(α) is the relative rotating velocity due to the presence of interactions (for the detailed calculation of ω(α) and expression of (27) the reader can refer to [19], [36] and [37]).

Results and discussion

The numerical computations were carried out under oscillatory shear flows (SAOS and LAOS) with a time step for the fiber position and orientation updates. Most simulations consider ∆t=0.005s. For all case studies a population of N=512 fibers (prolate spheroids) were distributed almost uniformly and isotropically within the REV. Using less than fibers does not allow a representative description of fiber population [41]. Their isotropic state was verified by applying equations (7), thus initially one gets:

a 2 1 3 0 0 0 1 3 0 0 0 1 3                             (28)    MathType@MTEF@5@5@+= feaahqart1ev3aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqiVu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaaeaaaaaaaaa8 qacaWGHbWdamaaBaaaleaapeGaaGOmaaWdaeqaaOWdbiabloKi7maa dmaapaqaauaabaqadmaaaeaapeWaaSaaa8aabaWdbiaaigdaa8aaba Wdbiaaiodaaaaapaqaa8qacaaIWaaapaqaa8qacaaIWaaapaqaa8qa caaIWaaapaqaa8qadaWcaaWdaeaapeGaaGymaaWdaeaapeGaaG4maa aaa8aabaWdbiaaicdaa8aabaWdbiaaicdaa8aabaWdbiaaicdaa8aa baWdbmaalaaapaqaa8qacaaIXaaapaqaa8qacaaIZaaaaaaaaiaawU facaGLDbaacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqG GaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabc cacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeii aiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeikaiaabkdacaqG4a GaaeykaiaabccacaqGGaGaaeiiaaaa@5CAB@

Table 1: Parameters defining the different case studies.

Case study

Total number of fibers

Length distribution (L)

1

512

512 with L = 30 mm

2

512

512 with L = 3 mm

3

512

256 with L1 = 30 mm (long fiber population)
256 with L2 = 3 mm (short fiber population)

The average aspect ratio of the prolate spheroids representing the fibers (equation (1)), can be obtained from the data reported in Table 1. In all cases, the fiber diameter is 1 mm. In case studies 1 and 2, we have all the fibers having a length of 30 mm and 3mm respectively.

Then in case 3, half of the fibers (i.e. fibers) are considered long fibers of length L1= 30mm and the other half short fibers of length L2= 3mm. The choice of those parameters are consistent with experimental studies made on fiber reinforced composites [42]. Furthermore, case 3 was chosen to have a study as close as possible to real industrial composites where fiber breakage might occur [24].

For a better, solid definition of "small" and "large" oscillatory flows, one can refer to recent literature that defined LAOS for polymeric liquids as operating for γ0 >1 [43]. Therefore, for the case studies of Table 1 an amplitude of γ0=1 and an amplitude of γ0=5 was applied respectively for SAOS and LAOS. The frequency of the oscillation was chosen such that w in equation (8) is fixed to 1 rad/s for all three case studies.

All three case studies were performed in the dilute and concentrated regimes. The concentrations were selected with respect to equations (2) and (3). For the dilute regime a volume fraction of φ=0.1% was considered for all case studies. Regarding the concentrated regime a volume fraction of φ=35% was considered for all case studies.

Kinematics evolution for case study 1 (L=30mm)

For the first case study, Figure 5 shows the time evolution of the a11, component with time when a SAOS is applied to the fiber suspensions, with the shear rate according to equation (8). a11 is plotted versus time.

Figure 5: Time evolution of the a11 component (SAOS: γ0 = 1; w = 1)

In the dilute regime there are no interactions between the fibers, therefore there are no restrictions for the fibers rotation and each fiber can rotate independently from its neighboring one. Thus, the period of a11 in Figure 5 has been checked to be equal to the period of γ ̇ (if the period is denoted by T, then T≃6.28 s from Figure 5). Moreover, when the flow is reversed (i.e. when γ < 0, Figure 6), the microstruture recovers the almost isotropic state in the absence of interactions.

In the concentrated regime, the evolution of the kinematics from Figure 5 is pseudo-periodic. The periodicity was broken due to the interactions taking place between fibers (contact and lubrication forces), that restricted the motion of the fibers with the same period of the applied shear rate.

Figure 6: Time evolution of γ ̇ (SAOS: 00 = 1; w = 1)

This effect, however, is not considerable enough due to the relatively low number of interactions, NI, occurring between the fibers as depicted in Figure 7. N1 is low, relative to the total number of fibers ( fibers) since the interactions are unilateral. Thus, fibers that were in contact can detach instantaneously when reversing the flow. When the flow is suddenly reversed, the number of interactions will diminish because fibers will tend to detach (Figure 7). Moreover, the loss in the number of interactions results in less entangled fibers.

When a LAOS is applied, the fibers got more aligned in the direction of the velocity field (equation (9)) as the magnitude of the shear rate increased from to . The shear rate being greater here, it forces the fibers to get closer to an aligned configuration ( from Figure 8). Although it may seem that higher shear rates induce more frequent interactions but as fibers become more aligned, the interactions decrease. This is why, here, was found to decrease in the same time intervals as when the fibers got more aligned (Figure 9). Because of the tendency of fibers to be, most of the time, more oriented in the flow direction (wider peaks for than the peaks observed in Figure 5), there aren’t enough time for the fiber-fiber interactions to happen (This is illustrated by the narrow peaks for in Figure 9). This is why the rotation of the fibers are not perturbed by interactions when the flow is reversed, and remains almost periodic in the concentrated regime.

Figure 7: Time evolution of the number of interactions (SAOS: γ0 = 1; w = 1)

Figure 8: Time evolution of the a11 component (LAOS: γ0 = 5; w = 1)

Figure 9: Time evolution of the number of interactions (LAOS: γ0 = 5; w = 1)

Kinematics evolution for case study 2(L = 3mm)

The same analysis made in section 3.1 can be made here for a SAOS applied in the dilute regime. However, in the concentrated regime, since all of the fibers have a lower aspect ratio (i.e. shorter fibers of length ), the system will be considered less concentrated with respect to equation (3). Therefore, the fiber-fiber interactions will be less considerable [39] than in the concentrated regime of case study 1, since the aspect ratio decreased by a factor of (r=3 while r=30 in case study 1). With almost no interactions occurring between the fibers (Figure-10), the time orientation evolution in the dilute and concentrated regimes are approximately the same and superimposed (Figure-11). When a LAOS is applied, the fibers reach a less aligned orientation state that the one in Figure 8 (a11 0.73 from Figure-12). This is expected since fibers are much shorter here. Again was found to be negligible and its evolution was approximately the same as the one found in Figure -10.

Kinematics evolution for case study 3 ( L1=30 mm and L2=3 mm )

Under SAOS, the kinematics time evolution in the polydisperse (i.e. different lengths of fibers) case study 3 is represented in Figure 13. The time evolution of the orientation, is similar to one found in Figures-5 and Figure- 11. In fact, when interactions are not considered (dilute regime), the solution exhibits a periodic evolution since which is expected. Introduction of fiber-fiber interactions (concentrated regime) was found to mostly affect the rotation of long fibers, while the kinematics of short fibers seem unperturbed by the activation of fiber-fiber interactions. Indeed by having of fibers of length L1= 30mm it was found that interactions are only due to the existence of long fibers. is depicted in Figure 14, where it is in average less the number of interactions obtained in case study 1. This is logical, since only half of the population is characterized with a high aspect ratio of r=30.

Figure 10: Time evolution of the number of interactions (SAOS: γ0= 1; w = 1)

Under LAOS, the fibers will be more aligned in the direction of the velocity field. But, the short fiber population here ( L2=3mm) with smaller aspect ratio produces sporadic rotations. As a result, a peak of a11=0.78 is reported (Figure-15) which is less than what was evidenced in Figure 8. Consequently, the number of interactions here, during the the time intervals where the fibers were mostly aligned in the direction of the flow, was found to be greater than the one found in Figure-9.

Figure 11: Time evolution of the a11 component (SAOS: γ0 = 1; w = 1)

Figure 12: Time evolution of the a11 component (LAOS: γ0 = 5; w = 1)

Figure 13: Time evolution of the a11 component (SAOS: γ0 = 1; w = 1)

Conclusions

Being used to investigate polymeric viscoelastic properties, oscillatory shear flows are known frequently to deliver unpredictable fiber orientations experimentally. In this paper, direct numerical simulations were proposed to investigate the time evolution of the orientation and kinematics of fiber suspensions under small and large amplitude oscillatory shear flow. In the dilute regime (under SAOS and LAOS), and for all cases considered in Table 1, the kinematic time evolution of the fibers (in terms of the component of the orientation tensor) was found to be periodic as the fibers, free of interactions, were rotating and oscillating in phase with the same period of the applied shear rate.

In the concentrated regime (under SAOS and LAOS) it was found for all case studies, that relative to the number of fibers, the number of interactions was low due to the fibers that kept detaching when the flow was reversed. When the SAOS was applied to cases and of Table 1 the periodicity in the time evolution of the component was slightly perturbed, due to those interactions. In case it was found that there were approximately no interactions under SAOS and LAOS, given the fact that the fibers aspect ratio decreased by a factor of .

Figure 14: Time evolution of the number of interactions (SAOS: γ0 = 1; w = 1)

Figure 15: Time evolution of the a11 component (LAOS: γ0 = 5; w = 1)

Figure 16: Time evolution of the number of interactions (LAOS: γ0 = 5; w = 1)

When the LAOS was applied to cases , and (concentrated regime) of Table 1 the time evolution of the component remained almost periodic and was not effected by fiber-fiber interactions. This was due, to fibers orienting more in the direction of the velocity field, and becoming more parallel to each others, thus decreasing the chance of them being in interaction.

The objective of having a fiber population where long fibers and short fibers are mixed together, was to try to mimic real industrial composites processes where fiber breakage might occur [24], therefore restricting the alignment of the fibers in a preferential direction. This was evidenced in case , where the short fibers restricted the component of the orientation tensor to no more than 0.78 .

The perspectives of this work is to evaluate the rheological macroscopic properties, such as stress, storage/loss moduli, and viscosity, for the fiber suspensions under SAOS and LAOS; this would be possible since those properties depend directly on the investigated second order orientation tensor, herein. This would be interesting to explore, since studies show that fiber suspensions in a polymer matrix manifest remarkable non-linear viscoelasticity even at very low deformation amplitude [44], in the framework of concentrated regimes. Furthermore, increasing the fiber content, can increase the storage and loss moduli of the neat matrix/resin [44]. Moreover, this will allow one to compare the macroscopic properties computed through DNS simulations, with the ones reported through experimental measurements.

References

  1. HPS  Abdul  Khalil,  AH  Bhat,  A  Abu  Bakar,  Paridah  Md  Tahir,  ISM  Zaidul,  and   M Jawaid. Cellulosic nanocomposites from natural fibers for medical applications: A review. In Handbook of polymer nanocomposites. Processing, performance and application.2015;475–511.
  2. Advani SG. Flow and rheology in polymer composites manufacturing, volume 10. Elsevier science. 1994.
  3. Advani SG, E Murat Sozer, and L Mishnaevsky Jr.  Process modeling in composites manufacturing. Appl. Mech. Rev. 2003;56(5):B69–B70.
  4. Tucker CL and Dessenberger RB. Flow and rheology in polymer composites manufacturing. Flow and rheology in polymer composites manufacturing. 1994.
  5. Randy S Bay and Charles L Tucker  III.  Fiber orientation in simple injection moldings. part i: Theory and numerical methods. Polymer composites.1992;13(4):317–331.
  6. Du-Hwan Chung and Tai-Hun Kwon. Fiber orientation in the processing of polymer composites. Korea-Australia Rheology Journal.2002;14(4):175–188.
  7. Steven M Dinh and Robert C Armstrong. A rheological equation of state for semicon- centrated fiber suspensions. Journal of Rheology. 1984;28(3):207–227.
  8. GG Lipscomb II, Morton M Denn, DU Hur, and David V Boger. The flow of fiber suspen- sions in complex geometries. Journal of Non-Newtonian Fluid Mechanics. 1988;26(3):297– 325.
  9. Christopher JS Petrie. The rheology of fibre suspensions. Journal of Non-Newtonian Fluid Mechanics.1999;87(2-3):369–402.
  10. Sepehr M, Ausias G and Carreau PJ.  Rheological properties  of  short  fiber filled polypropylene in transient shear flow. Journal of Non-Newtonian Fluid Mechanics. 2004;123(1):19–32.
  11. MandDoi and SF Edwards.  Dynamics of rod-like  macromolecules  in  concentrated  solution. part 1. Journal of the Chemical Society, Faraday Transactions 2: Molecular and Chemical Physics.1978;74:560–570.
  12. George Barker Jeffery. The motion of ellipsoidal particles immersed in a viscous fluid.Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character.1922;102(715):161–179.
  13. Abd El-Rahman, I Ahmed, and Charles L Tucker. Mechanics of random discontinuous long-fiber thermoplastics—part i: Generation and characterization of initial geometry. Journal of Applied Mechanics. 2013;80(5):051007.
  14. AI Abd El-Rahman and CL Tucker III. Mechanics of random discontinuous long-fiber thermoplastics. Part II: Direct simulation of uniaxial compression. Journal of Rheology. 2013;57(5):1463–1489.
  15. Folgar F and Tucker CL III. Orientation behavior of fibers in concentrated suspensions. Journal of reinforced plastics and composites. 1984;3(2):98–119.
  16. Batchelor GK.  The stress generated in a non-dilute suspension of elongated particles by pure straining motion. Journal of Fluid Mechanics.1971;46(4):813–829.
  17. Toll S and Jan-Anders E M.  Dynamics of a planar concentrated fiber sus- pension with non-hydrodynamic interaction. Journal of Rheology. 1994;38(4):985–997.
  18. Le Corre S, Dumont P, Orgeas L and Favier D. Rheology of highly concentrated planar fiber suspensions. Journal of Rheology. 2005;49(5):1029–1058.
  19. Ausias G, Fan XJ, and Tanner RI. Direct simulation for concentrated fibre suspensions in transient and steady state shear flows. Journal of non-newtonian fluid mechanics.2006;135(1):46–57.
  20. John D Ferry. Viscoelastic properties of polymers. John Wiley & Sons, 1980.
  21. Kitano T, Kataoka T, and Nagatsuka Y. Dynamic flow properties of vinylonfibre and glass fiber reinforced polyethylene melts. Rheologicaacta. 1984;23(4):408–416.
  22. AT Mutel and MR Kamal. Characterization of the rheological behavior of  fiberfilled polypropylene melts under steady and oscillatory shear using cone-and-plate and rotational parallel platerheometry. Polymer composites. 1986;7(5):283–294.
  23. Joseph P Greene and James O Wilkes. Steady-state and dynamic properties of concentrated fiber-filled thermoplastics. Polymer Engineering & Science. 1995;35(21):1670–1681.
  24. Mobuchon C, Carreau PJ, Heuzey MC, Sepehr M, and Ausias G. Shear and extensional properties of short glass fiber reinforced polypropylene. Polymer composites.2005;26(3):247–264.
  25. Corre SL, Caillerie D, Orgeas L, and Favier D.  Behavior of a net of fibers linked by viscous interactions: theory and mechanical properties.  Journal of the Mechanics and Physics of Solids.2004;52(2):395–421.
  26. JJ  Dumont P,  Corre SL,  Orgeas L,  and  Favier D.   A numerical analysis of the evolution of bundle orientation in concentrated fibrebundle suspensions. Journal of non-newtonian fluid mechanics. 20009;160(2-3):76–92.
  27. Yamamoto S and Matsuoka T. Dynamic simulation of a platelike particle dispersed system. The Journal of chemical physics. 1997;107(8):3300–3308.
  28. C Pozrikidis. Orientation statistics and effective viscosity of suspensions of elongated particles in simple shear flow. European Journal of Mechanics-B/Fluids. 2005;24(2):125–136.
  29. Sundararajakumar RR and Koch DL. Structure and properties of sheared fiber suspensions with mechanical contacts. Journal of Non-Newtonian Fluid Mechanics.1997;73(3):205–239.
  30. Satoru Yamamoto and Takaaki Matsuoka. A method for dynamic simulation of rigid and flexible fibers in a flow field. The Journal of chemical physics.1993;98(1):644–650.
  31. Advani SGand Charles L Tucker III.  The use of tensors to describe and predict fiber orientation in short fiber composites. Journal of rheology. 1987;31(8):751–784.
  32. Bossis G, Meunier A, and Sherwood JD. Stokesian dynamics simulations of particle trajectories near a plane. Physics of Fluids A: Fluid Dynamics. 1991;3(8):1853–1858.
  33. Born M and Karman V. Overshoots in space lattices. J Phys.1912;13:297–309.
  34. Allen MP, Tildesley DJ. Computer simulation of liquids. 1987.
  35. Townsend AK and Wilson HJ.  Small-and large-amplitude oscillatory rheome- try with bead–spring dumbbells in stokesian dynamics to mimic viscoelasticity.  Journal   of Non-Newtonian Fluid Mechanics.2018;261:136–152.
  36. S. Kim and S. J. Karrila. Microdynamics: Principles and Selected Applications. Butterworth–Heinemann, Boston, 1991.
  37. RabihMezher. Modeling and Simulation of concentrated suspensions of short, rigid and flexible fibers. PhD thesis. 2015.
  38. Charles L Tucker III. Flow regimes for fiber suspensions in narrow gaps. Journal of Non-Newtonian fluid mechanics. 1991;39(3):239–268.
  39. Mezher R, Abisset-Chavanne E, Ferec J, Ausias G, and Chinesta F. Direct simulation of concentrated fiber suspensions subjected to bending effects. Modelling and Simulation in Materials Science and Engineering.2015;23(5):055007.
  40. Harlen OG, Sundararajakumar RR and Koch DL.  Numerical simulations of a sphere settling through a suspension of neutrally buoyant fibres. Journal of Fluid Mechanics.1999;388:355–388.
  41. Chavanne EA, Mezher R, Corre SL, Ammar A and Chinesta F. Kinetic theory microstructure modeling in concentrated suspensions. Entropy.2013;15(7):2805–2832.
  42. Tassew STand Lubell AS. Mechanical properties of glass fiber reinforced ceramic concrete. Construction and Building Materials.2014;51:215–224.
  43. Giacomin AJ, Bird RB, Johnson LM and Mix AW.  Large-amplitude oscillatory shear flow from the corotationalmaxwell model. Journal of Non-Newtonian Fluid Mechanics.2011;166(19-20):1081–1099.
  44. Ferec J, Heuzey MC, Ausias G, and Carreau PJ.  Rheological behavior of fiber-filled polymers under large amplitude oscillatory shear flow. Journal of non-newtonian fluid mechanics. 2008;151(1-3):89–100.
//