# 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:

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:

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

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:

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}^{\left(\alpha \right)}={p}_{1}^{\left(\alpha \right)}x+{p}_{2}^{\left(\alpha \right)}y+{p}_{3}^{\left(\alpha \right)}z$ 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].

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:

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 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]:

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:

(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:

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:

The rate of strain tensor is written as:

and the vorticity tensor is given by:

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(α) 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:

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:

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:

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]:

### Interaction forces

The resulting interaction force on fiber α is:

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]:

with:

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]:

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:

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:

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 . The two new resistance tensors ξ(α) and ${\mathcal{H}}^{\left(\alpha \right)}$ , 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:

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:

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

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:

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.