During processing of fiber composites, the fiber-induced stresses influence the local flow fields, which, in turn, influence the stress distribution and the fiber orientation. Therefore, it is crucial to be able to predict the rheology of fiber-filled polymer composites. In this study, we investigate the fiber orientation kinetics and rheological properties of fiber composites in uniaxial extensional flow by comparing direct numerical finite element simulations to experimental results from our previous study [Egelmeers *et al.*, “In-situ experimental investigation of fiber orientation kinetics during uniaxial extensional flow of polymer composites,” J. Rheol. **68**, 171–185 (2023)]. In the simulations, fiber–fiber interactions only occur hydrodynamically and lubrication stresses are fully resolved by using adaptive meshing. We employed a 7-mode and a 5-mode viscoelastic Giesekus material model to describe the behavior of, respectively, a strain hardening low-density polyethylene (LDPE) matrix and a non-strain hardening linear LDPE matrix, and investigated the influence of the Weissenberg number, strain hardening, and fiber volume fraction on the fiber orientation kinetics. We found that none of these parameters influence the fiber orientation kinetics, which agrees with our experimental data. The transient uniaxial extensional viscosity of a fiber-filled polymer suspension is investigated by comparing finite element simulations to a constitutive model proposed by Hinch and Leal [“Time-dependent shear flows of a suspension of particles with weak Brownian rotations,” J. Fluid Mech. **57**(4), 753–767 (1973)] and to experimental results obtained in our previous study [Egelmeers *et al.*, “In-situ experimental investigation of fiber orientation kinetics during uniaxial extensional flow of polymer composites,” J. Rheol. **68**, 171–185 (2023)]. The simulations describe the experimental data well. Moreover, high agreement is found for the transient viscosity as a function of fiber orientation between the model and the simulations. At high strains for high fiber volume fractions, however, the simulations show additional strain hardening, which we attribute to local changes in microstructure.

## I. INTRODUCTION

Fiber-filled polymer composites are important materials for the high tech industry because they combine important material properties, e.g., high strength, low weight, and low cost of manufacturing. However, fiber-filled polymers offer distinct challenges: the fiber orientation induced during processing influences the mechanical properties, and the anisotropy of the fibers can cause warpage during the cooling step in the manufacturing process of a product. Modeling the fiber orientation kinetics and the rheological properties during processing is key to mitigating these limitations.

Modeling of fiber orientation can be performed on different length scales: the micro scale,^{3–9} where individual fibers are considered, and the macro, or continuum scale,^{10–15} which is related to the length scale of an actual fiber-filled part. At the macro scale, fiber orientation is commonly modeled using the Folgar–Tucker model or variations thereof.^{10–12} The basis of these models is Jeffery's equation,^{16} which describes the orientation kinetics of a single ellipsoidal particle suspended in a Newtonian fluid. Other terms in the Folgar–Tucker models are phenomenological and describe, e.g., rotary diffusion, which is related to Brownian motion and slows down the fiber orientation. In practice, the parameters in the additional terms of the Folgar–Tucker models are determined by fitting experimental data to include effects of interparticle interactions.^{17} The Folgar–Tucker models are computationally tractable, but the coarse-graining leads to a loss of information at the micro-scale.

At the micro-scale, different modeling approaches exist, e.g., particle-based methods,^{3–6} where the polymer matrix and the fibers are treated as particles, and grid-based methods,^{7–9} where the polymer matrix is treated as a continuous medium and the fibers as particles. Within these methods, a distinction can be made between the manners in which fluid–fiber interactions are treated, which can be either one-way coupled or two-way coupled. One-way coupling of the fluid–fiber interactions only considers the hydrodynamic forces exerted by the fluid on the fiber, and the fluid velocity is undisturbed by the fiber motion. Two-way coupling includes the backcoupling of the fiber motion on the fluid velocity and, in this way, also includes hydrodynamic long-range and short-range fiber–fiber interaction forces.^{18}

As Kruger *et al.*^{18} pointed out, existing micro-scale simulations mostly focus on shear flow. A limited number of investigations are performed for extensional flows. El-Rahman and Tucker^{19} considered fiber orientation in uniaxial compression. They neglected the presence of the matrix fluid and modeled the fibers as 3 D Timoshenko beam elements, and concluded that the overall fiber orientation is affine. Yamanoi and Maia^{5} investigated the fiber orientation kinetics in a Newtonian fluid for uniaxial extensional flow. They performed Stokesian dynamics simulations, while adding artificial fiber–fiber interaction forces to keep fibers from colliding, and investigated the influence of fiber aspect ratio, fiber bending flexibility, and fiber volume fraction on the fiber orientation kinetics. They concluded that these parameters do not influence the orientation tensor component in the flow direction. However, it remains unclear how viscoelastic matrix properties influence the fiber orientation in uniaxial extension.

The fiber orientation strongly depends on the flow fields applied during processing. The flow fields are in turn a function of the rheology.^{20,21} Therefore, understanding the rheology of fiber-filled polymer suspensions is of great importance. A significant amount of research has already been conducted to understand the rheology of suspensions with spherical particles, as comprehensively reviewed by Tanner.^{22,23} For fiber-filled suspensions, the fiber aspect ratio introduces additional complexity and challenges. Current theoretical models describing the rheology of fiber-filled suspensions often use slender-body theory.^{24–26} Some of these models do not take fiber orientation kinetics into account and assume a constant fiber orientation state.^{24,25} An example is an expression developed by Batchelor,^{24} who considered the semi-concentrated regime and assumed that all fibers are aligned in the stretching direction. His expression is based on a parameter that characterizes the distance from a fiber to its nearest neighbor, which depends on the fiber volume fraction. Similar expressions are derived for the dilute and concentrated regime by Petrie^{25} and Pipes *et al.*,^{27,28} respectively. Other rheological models include fiber orientation kinetics.^{26,29–31} In these models, the Folgar–Tucker equation is solved to compute the fiber orientation as a function of time and the composite's uniaxial extensional viscosity is a function of the fiber orientation. For fibers aligned in the stretch direction, the extensional viscosity is larger compared to that for a random fiber orientation. In uniaxial extension, these models will predict orientation-induced strain hardening. Yet, experimental studies investigating the rheology of fiber suspensions performed by Egelmeers *et al.*,^{1} Wagner *et al.*,^{32} and Wang *et al.*^{33} did not show the orientation-induced strain hardening. In contrast with these findings, Ferec *et al.*^{30} investigated the rheological properties of a glass fiber-filled polypropylene and did observe orientation-induced strain hardening.

As most rheological models assume a slender-body theory, a deeper understanding about the rheological properties of fiber-filled suspensions can be provided by direct numerical simulations. Stokesian dynamics simulations of fibers in a Newtonian matrix for uniaxial extension were performed by Yamanoi and Maia^{5} and also showed strain hardening. However, it still remains unclear how the complex interplay between a viscoelastic matrix, fiber orientation, and fiber–fiber interactions influences the rheological properties in a fiber-filled polymer in uniaxial extensional flow.

In this study, direct numerical simulations for uniaxial extensional flows of fiber suspensions with a viscoelastic matrix will be performed using the finite element method (FEM). The fiber orientation kinetics and the corresponding composites' rheology are studied as a function of the Weissenberg number, strain hardening, and fiber volume fraction, and will be compared to experimental data obtained in a previous study by Egelmeers *et al.*^{1} The rheology resulting from the FEM simulations will be compared to the rheology from a constitutive model proposed by Hinch and Leal.^{2} The FEM simulations include the backcoupling from the fiber motion on the fluid motion and consider the case where the interaction between particles is purely hydrodynamical in nature and direct particle–particle interaction forces (e.g., electrostatic or van der Waals) are negligible. Moreover, by fully resolving the lubrication stresses in between particles, fiber collision is prevented without the need for additional numerical measures to prevent overlap.

## II. PROBLEM DESCRIPTION

**) is equal to $ u * + u \u0302$, where $ u *$ is the periodic part of the velocity, and $ u \u0302$ the velocities described by macroscopic uniaxial extensional field, according to**

*u**x*,

*y*, and

*z*are defined in Fig. 1. Note, in this manuscript, the notation $ [ . . ] T$ is used to indicate the Cartesian vector components. In the following, a triperiodic domain is identified by the integers

*i*,

*j*, and

*k*with the reference domain being

*i*= 0,

*j*= 0, and

*k*= 0. Between a point in the reference domain and its periodic image(s) in the other domains, the velocity difference equals

^{9,34}

*L*,

_{x}*L*, and

_{y}*L*will be used to indicate the time-dependent domain sizes.

_{z}### A. Particles in the triperiodic domain

*i*th particle. The angle

*θ*is called the alignment angle, which describes the angle between the particle orientation vector and the stretch direction, in this case the

*z*axis. If

*θ*is close to zero the particle is aligned with the

*z*axis. The second angle,

*ψ,*is the angle between the

*y*axis and the projection of the particle orientation vector on the

*xy*plane. Figure 2(a) visualizes a triperiodic rectangular box with two complete particles. One complete particle is “cut” into four parts by the domain boundaries. Figure 2(b) visualizes the consistency of the triperiodic domain and shows that the cut particles indeed make up one complete particle. The velocities of the individual cut particles can be corrected using Eq. (2).

### B. Balance equations

**is the fluid velocity vector and $\sigma $ is the Cauchy stress tensor**

*u**p*is the pressure,

**the unit tensor, $ \eta sol$ the Newtonian solvent viscosity, and $ D = ( \u2207 u + ( \u2207 u ) T ) / 2$ represents the rate of deformation tensor. The term $ 2 \eta sol D$ is the Newtonian stress tensor and $\tau $ is the viscoelastic stress tensor. The viscoelastic stress tensor for the Giesekus constitutive model in conformation tensor representation is given by**

*I*^{35}

**is the conformation tensor, $ \eta p$ the polymer viscosity, and**

*c**λ*is the polymer relaxation time. The latter two are related through the polymer modulus $ G = \eta p / \lambda $. Furthermore,

*α*is a (dimensionless) constitutive parameter, which influences shear thinning behavior and the symbol ( $ \u25bd$) denotes the upper-convected time derivative

### C. Boundary and initial conditions

**the outward pointed normal on surface $ \Gamma i$, continuity of the stress yields**

*n**et al.*

^{9}In this work, the latter approach has been chosen.

*i*th particle surface, $ n i$ the outwardly directed unit normal vector on that particle surface,

**are coordinates on that particle surface, and $ X i = [ X i , Y i , Z i ] T$ is the local center point coordinate of particle**

*x**i*. Assuming no-slip, the boundary condition for the velocity on the complete particle boundary $ \u2202 F i$ will therefore become

*i*th particle, and $ U i = [ U x , i , U y , i , U z , i ] T$ is the unknown local center point velocity of particle

*i*, except for the first particle where $ U 1 = 0$ is prescribed. For cut particles, different parts of the surface are associated with different periodic center points, which need to be taken into account in Eqs. (18) and (19). Moreover, in this case, the velocity difference for different periodic center points, $ U i$, is corrected using Eq. (2).

### D. Arbitrary Lagrange–Euler (ALE) formulation

^{37}The movement of the mesh has to be compensated, which only applies to the evolution equation of the conformation tensor, since this is the only equation that contains a convective term, which is rewritten as

### E. Bulk stresses

^{38}

^{,}

*V*is the volume of the domain, Ω is the volume occupied by the fluid in the domain, and $ F = F 1 \u222a F 2 \u222a \u2026 F N$ is the volume occupied by the particles in the domain, see Fig. 2(a). By substituting $ \sigma = \u2207 \xb7 ( \sigma T x )$ and afterward using the divergence theorem, the integral over the particles' volume can be rewritten to an integral over the particles' boundary surface, indicated by $ \u2202 F$. Furthermore, Eq. (7) is substituted in Eq. (22),

^{9,39}and finally, for a correct implementation of the cut particles' contribution to the volume-averaged stress, the local midpoint of the cut particle is subtracted from the surface coordinate of the particles, according to

*N*is the number of particles in the domain. Finally, the volume-averaged viscosity is computed, according to

### F. Anisotropy

*et al.*

^{1}The average particle orientation is computed using the particle orientation tensor

^{10}

*z*axis with initial random particle orientation, all off-diagonal components of

**at every time step are almost zero. Therefore, the diagonal components of**

*P***equal the eigenvalues of the tensor. The anisotropy is determined by comparing the average orientation components of this tensor in the**

*P**x*and

*y*directions with the average orientation component in the

*z*direction

*ν*are the eigenvalues of the particle orientation tensor. For a random particle orientation, $ P x x = P y y = P z z = 1 / 3$, so that the anisotropy, $ \xi s = \u2212 0.41$ and when all particles are perfectly aligned in the stretch direction, $ P x x = P y y = 0$ and

_{z}*P*= 1, so that $ \xi s = 1$.

_{zz}### G. RVE statistics

^{9}

*i*th realization. The interval that contains the true average with 95% confidence equals $ \xi \xaf s \u2009\xb1\u2009 2 s \xi \xaf s$ and $ \eta \xaf E + \u2009\xb1\u2009 2 s \eta \xaf E +$. In this study, it is ensured that after a strain of 2, the width of this interval does not exceed 8% of the average. Therefore, in the remainder of the paper, the number of realizations per simulation setting is such that the following criteria are met:

## III. NUMERICAL METHOD

The balance equations are solved using an in-house developed software package utilizing the finite element method (FEM). Discrete Elastic-Viscous Stress Splitting using the G approach (DEVSS-G), Streamline Upward Petrov Gaelerkin (SUPG), and the log representation of the conformation tensor are used to increase the stability of the simulations.^{40–44}

### A. Weak formulation

*q*, $ v$, $ V i , \u2009 \mu i$

**, $ \chi i$,**

*d***. Moreover, $ G = \u2207 u T$ is the velocity gradient, $ D v = 1 2 ( \u2207 v + ( \u2207 v ) T ) , \u2009 ( \xb7 , \xb7 )$ denotes the inner product in domain Ω, $ \u27e8 \xb7 , \xb7 \u27e9 \u2202 F$ denotes the inner product on particle surface $ \u2202 F$,**

*H**τ*and

*ν*are parameters for the SUPG and DEVSS-G stabilization, $ s = log \u2009 c$, and $ h ( G , s )$ is a tensor function.

^{43}

### B. Spatial discretization

Isoparametric tetrahedral Taylor–Hood elements are used for the velocity and the pressure, employing quadratic interpolation of the velocities (*P*_{2}) and linear interpolation for the pressures (*P*_{1}). For the velocity gradient and the conformation tensor, *P*_{1} interpolation is used. The mesh is generated using Gmsh.^{45} Adaptive local mesh refinement is implemented in between approaching particles. The linear system resulting from the set of equations outlined above is solved using Pardiso,^{46–48} which is a parallel direct solver integrated in the Math Kernel Library (MKL).

### C. Time discretization

^{49}In addition, to not be limited in time step size by very fast modes, the terms $ s n + 1 / \lambda $ and $ s \u0302 n + 1 / \lambda $ are added on the left and right sides of the equation, respectively. The second-order semi-implicit Gear scheme in log form now becomes

*n*and the value of a variable at the current time step with a superscript notation: i.e., $ u | t = n \Delta t = u n$. Moreover, $ s \u0302 n + 1 = 2 s n \u2212 s n \u2212 1$ is the second-order predictor of the log of the conformation tensor at the new time step. For relaxation times much smaller than the time step, the time-dependent terms will become small compared to the relaxation terms. The corresponding mode will therefore “step over the transient” into a steady state, with the local steady-state stress depending on the local Weissenberg number, $ \lambda | D |$.

^{50}

^{,}

*W*a matrix given by

*i*denotes the particle number. Initially, the quaternions for all particles are initialized as $ q | t = 0 = q 0 = [ 1 , 0 , 0 , 0 ] T$ and they are updated in the first step using an Euler integration scheme

*n*, to the particle orientation at the new time step,

*n*+ 1, for the

*i*th particle, which is done using a two-step procedure: First, the quaternion at the current time step is mapped to the initial state (using orthogonality of

*R*, i.e., $ R \u2212 1 = R T$), and subsequently, the initial state is mapped to the new time step, according to

^{51}

**. The orientation of the particle is updated according to**

*q***is a diagonal tensor with $ A x x = A y y = exp \u2009 ( \u2212 \epsilon \u0307 H \Delta t / 2 )$ and $ A z z = exp \u2009 ( \epsilon \u0307 H \Delta t )$. Note that nodes shared by particles and periodic surfaces are updated according to the particle displacement.**

*A*^{37}

*k*is a diffusion coefficient proportional to the inverse of the element size allowing the coarser elements away from the surfaces to deform, while the finer elements close to the surfaces are translated, to minimize element deformation.

_{e}### D. Remeshing and projection

^{37}

Meshing is, as explained in Sec. III B, performed using Gmsh. After meshing, a projection problem^{52} is solved to determine the solution of the new mesh in the previous time step ( $ c n , \u2009 c n \u2212 1 , \u2009 u n , \u2009 u n \u2212 1$). The ALE formulation also requires to find the coordinates of the new mesh at the previous time step.

In case Gmsh is not able to successfully generate a mesh, the vertices of the triperiodic domain are translated, while the locations of the particles are maintained. Then, the meshing and projection procedure as explained above are repeated. However, the old and the new domains do not necessarily overlap, and this is solved by finding the projection solution in the corresponding periodic image, taking the shift in periodic velocities Eq. (2) into account. Figure 3 shows the schematic remeshing and projection procedure in the case of translated domain vertices.

## IV. NUMERICAL MODEL VERIFICATION

In this section, different studies are performed to verify the numerical implementation. First, mesh and time convergences tests are performed for a cut particle in a Newtonian matrix. Afterward, mesh convergence is investigated for particles in a viscoelastic matrix.

### A. Convergence study

^{−1}is applied in uniaxial extensional flow. The initial orientation vector of the particle equals $ p = [ 1 2 2 , 0 , 1 2 2 ] T$. The particle and the domain are chosen such that the particle is cut into four parts, see Fig. 4(a). The domain will be extended in

*z*direction. The finite element simulations are compared to the analytical equation describing the orientation of a single spheroid in a Newtonian fluid for uniaxial extensional flow.

^{16,53}For a particle with initial orientation in the

*xz*plane, the analytical expressions for the rotational velocity and alignment angle are, respectively,

*L*/

*D*). The alignment angle

*θ*describes the angle between the stretch direction and the fiber orientation vector ( $ \theta = arccos ( p z )$). The evolution of the other angle to fully describe the orientation of the particle,

*ψ*, is not of interest since $ \psi ( t ) = 0$.

*y*-component of the rotational velocity is computed in the FEM simulation using the boundary conditions as described above. The relative error is defined as

*L*is the length of the prolate spheroid (

*L*= 4). In the triperiodic FEM simulations, the particle–particle distance equals the domain size, which introduces an error compared to the analytical expression that assumes no particle–particle interactions. Therefore, different domain sizes are considered in the mesh convergence plot. For the smallest domain size ( $ L 0 , x / L = 2$), this error contribution is clearly visible; for the other domain sizes, the error contribution due to particle–particle interactions becomes dominant at smaller element sizes. Moreover, in the region where the error is dominated by the mesh size, a slope of −3 is visible, which is expected for the second-order Taylor–Hood elements used here.

### B. Determination of the mesh size

In this section, viscoelastic particle orientation simulations are performed for spherocylindrical particles. The goal of this section is to determine the mesh size characteristics at which the anisotropy of the fibers in the bulk and the transient uniaxial extensional viscosity are converged. A single-mode Giesekus material model, with $ \eta p = 10 5$ Pa s, *λ* = 1 s, $ \alpha = 0.005$, and $ \eta sol = 0$ Pa s, is used in these simulations. Furthermore, a Hencky strain rate of 1 s^{−1} is applied with a time step of 0.01 s until a linear strain of 4 is reached. Every simulation in this section uses the same initial particle orientation and initial particle positions. The RVE consists of 16 particles with aspect ratio 10 and the particle volume fraction is 10%. Five different mesh size characteristics are considered. There are two parameters that control the mesh size: the number of nodes on the diameter of the particle, $ n d$, and the number of nodes between particles, $ n f \u2212 f$. Table I indicates the mesh size characteristics used for each mesh, and it includes the number of elements in the mesh initially and at the end of the simulation.

Mesh ID . | $ n d$ . | $ n f \u2212 f$ . | $ n elem , start$ . | $ n elem , end$ . |
---|---|---|---|---|

M1 | 6 | 1 | 22.000 | 47.000 |

M2 | 8 | 1 | 35.000 | 57.500 |

M3 | 10 | 2 | 54.000 | 110.000 |

M4 | 12 | 2 | 75.000 | 135.000 |

M5 | 12 | 3 | 80.000 | 190.000 |

Mesh ID . | $ n d$ . | $ n f \u2212 f$ . | $ n elem , start$ . | $ n elem , end$ . |
---|---|---|---|---|

M1 | 6 | 1 | 22.000 | 47.000 |

M2 | 8 | 1 | 35.000 | 57.500 |

M3 | 10 | 2 | 54.000 | 110.000 |

M4 | 12 | 2 | 75.000 | 135.000 |

M5 | 12 | 3 | 80.000 | 190.000 |

The different meshes are shown in Fig. 5. Every image shows a clipping plane in the mesh. Comparing the particles itself in Fig. 5(a) with the particles in Fig. 5(e), one clearly sees the decrease in element size on the particles' surface mesh for increasing $ n d$. Also, the influence of mesh size parameter, $ n f \u2212 f$, is clearly visible by comparing the area between the two particles at the right of Figs. 5(a) and 5(e).

Simulations are performed for these different meshes to compute the anisotropy of the fibers in the bulk and the transient uniaxial extensional viscosity. The results are shown in Fig. 6. The figure clearly shows that for the considered range of mesh size characteristics, there is no significant influence of mesh size characteristics on both bulk properties, as all lines coincide. To show that there are local differences for the different meshes, the trace of the conformation tensor, computed at a cross section of the RVE, is plotted for the different meshes at a linear strain of 4 in Fig. 7. It can be observed that the trace of the conformation tensor increases at the ends of the particles and when particles are located exactly underneath each other. For the remaining simulations in this study, mesh size characteristic M3 is chosen, which optimizes robust simulations, accuracy and computational efficiency. Moreover, the anisotropy and the bulk stress are fully converged as shown in Fig. 6.

## V. CONSTITUTIVE MODEL FOR RHEOLOGY OF FIBER SUSPENSIONS

^{2}

^{,}

^{17}For a completely random initial fiber orientation, $ P 0$ is a diagonal tensor with $ P x x = P y y = P z z = 1 / 3$, a natural closure approximation yields the most accurate predictions, whereas for a pseudorandom initial fiber orientation, the quadratic closure approximation ( $ P = PP$) is more accurate, which is demonstrated in Sec. S1 of the supplementary material. The fiber orientation tensor

**is a function of time. Its evolution is computed using Jeffery's equation in tensor form**

*P*^{17}

*β*is the particle shape factor, which is related to geometry and aspect ratio. For the particles in this work,

*β*= 1. By combining Eqs. (55) and (56), predictions can be made for the transient uniaxial extensional viscosity. These predictions will be compared to the predictions of the FEM simulations, where the viscosity can be obtained numerically (as explained in Sec. II E).

## VI. RESULTS

In this section, the fiber orientation kinetics and the transient uniaxial extensional viscosity are numerically investigated and the results are compared to experimental results obtained by Egelmeers *et al.*^{1} First, constitutive material model parameters are determined for the matrix to ensure consistency with the experiments.

### A. Matrix constitutive model parameters

In our previous study, the fiber orientation kinetics in extension were experimentally investigated.^{1} In this previous study, the matrix of the experimental system consisted either of a strain-hardening low-density polyethylene (LDPE), or a non-strain-hardening linear low-density polyethylene (LLDPE). In the current work, the rheology of both materials is described by a multi-mode Giesekus constitutive material model. The relaxation time and polymer viscosity for each mode were determined from small-amplitude oscillatory shear (SAOS) measurements in our previous work.^{1} In the current work, the remaining non-linear Giesekus model parameters *α _{i}* for each mode were determined using a fitting procedure on the transient shear viscosity and transient uniaxial extensional viscosity data. The resulting Giesekus model parameters that describe the rheology of the unfilled LDPE and LLDPE are summarized in Tables II and III, respectively.

Mode . | λ (s)
. | $ \eta p$ (Pa s) . | α
. |
---|---|---|---|

1 | $ 1.00 \xd7 10 \u2212 2$ | $ 9.29 \xd7 10 2$ | 0.49 |

2 | $ 8.51 \xd7 10 \u2212 2$ | $ 3.11 \xd7 10 3$ | 0.49 |

3 | $ 5.18 \xd7 10 \u2212 1$ | $ 8.91 \xd7 10 3$ | 0.49 |

4 | $ 3.02$ | $ 2.18 \xd7 10 4$ | 0.25 |

5 | $ 1.65 \xd7 10 1$ | $ 3.82 \xd7 10 4$ | 0.11 |

6 | $ 7.53 \xd7 10 1$ | $ 3.76 \xd7 10 4$ | 0.11 |

7 | $ 2.74 \xd7 10 2$ | $ 2.34 \xd7 10 4$ | 0.08 |

Mode . | λ (s)
. | $ \eta p$ (Pa s) . | α
. |
---|---|---|---|

1 | $ 1.00 \xd7 10 \u2212 2$ | $ 9.29 \xd7 10 2$ | 0.49 |

2 | $ 8.51 \xd7 10 \u2212 2$ | $ 3.11 \xd7 10 3$ | 0.49 |

3 | $ 5.18 \xd7 10 \u2212 1$ | $ 8.91 \xd7 10 3$ | 0.49 |

4 | $ 3.02$ | $ 2.18 \xd7 10 4$ | 0.25 |

5 | $ 1.65 \xd7 10 1$ | $ 3.82 \xd7 10 4$ | 0.11 |

6 | $ 7.53 \xd7 10 1$ | $ 3.76 \xd7 10 4$ | 0.11 |

7 | $ 2.74 \xd7 10 2$ | $ 2.34 \xd7 10 4$ | 0.08 |

Mode . | λ (s)
. | $ \eta p$ (Pa s) . | α
. |
---|---|---|---|

1 | $ 1.15 \u2009 \xd7 10 \u2212 2$ | $ 3.97 \xd7 10 3$ | 0.49 |

2 | $ 6.35 \u2009 \xd7 10 \u2212 2$ | $ 7.11 \xd7 10 3$ | 0.49 |

3 | $ 3.18 \u2009 \xd7 10 \u2212 1$ | $ 5.34 \xd7 10 3$ | 0.49 |

4 | $ 2.06$ | $ 1.61 \xd7 10 3$ | 0.49 |

5 | $ 1.70 \xd7 10 1$ | $ 5.44 \xd7 10 2$ | 0.49 |

Mode . | λ (s)
. | $ \eta p$ (Pa s) . | α
. |
---|---|---|---|

1 | $ 1.15 \u2009 \xd7 10 \u2212 2$ | $ 3.97 \xd7 10 3$ | 0.49 |

2 | $ 6.35 \u2009 \xd7 10 \u2212 2$ | $ 7.11 \xd7 10 3$ | 0.49 |

3 | $ 3.18 \u2009 \xd7 10 \u2212 1$ | $ 5.34 \xd7 10 3$ | 0.49 |

4 | $ 2.06$ | $ 1.61 \xd7 10 3$ | 0.49 |

5 | $ 1.70 \xd7 10 1$ | $ 5.44 \xd7 10 2$ | 0.49 |

Figures 8(a)–8(c) show the transient extensional viscosity, transient shear viscosity, and transient first normal stress difference, respectively, with experimental data obtained from our previous work (markers) and the Giesekus material model predictions (solid lines) with the parameters indicated in Table II for unfilled LDPE. The figures show that the Giesekus model prediction agrees well with the experimentally obtained transient behavior of unfilled LDPE in both extension and shear. These figures also show transient model predictions at rates of 10 and 100 s^{−1}, which were experimentally not measurable with the used setups, but relevant because high local rates occur between particles. To get experimental rheological information at these high rates, the Cox–Merz rule is applied on SAOS data obtained from our previous work.^{1} Fig. 8(d) shows the experimental complex viscosity data (markers), which, according to the Cox–Merz rule, are equivalent to steady shear data. The figure also includes the steady-state viscosity prediction as a function of the shear rate of the 7-mode Giesekus model describing LDPE (solid line). The model agrees reasonably well with the experimental data, showing that the multi-mode Giesekus model is able to predict the steady shear viscosity of LDPE well, also at high rates.

Figures 9(a)–9(c) show the transient extensional viscosity, transient shear viscosity, and transient first normal stress difference, respectively, with experimental data obtained in our previous work (markers) and the Giesekus material model prediction (solid lines) with the parameters indicated in Table III for unfilled LLDPE. The model agrees well with the transient behavior of LLDPE. Figure 9(d) shows the experimentally obtained complex viscosity data (markers), and the steady-state viscosity prediction of the 5-mode Giesekus model describing LLDPE (solid line). The model agrees well with the experimental data, implying that the multi-mode Giesekus model is able to predict the steady state shear viscosity of LLDPE, also at high rates.

### B. *In situ* deformation and stress fields

In this section, the deformation gradient around the fibers and the stresses in the matrix as a function of strain will be investigated. For this, a FEM simulation is performed with 16 fibers in the RVE, which have a fiber aspect ratio of 10. Furthermore, in this simulation, the fiber volume fraction is 10%, the applied Hencky strain rate is 1 s^{−1}, and a 7-mode Giesekus material model describing the behavior of LDPE, see Table II, is considered. All individual fibers for this simulation at different strains are shown in Fig. 10, where the color indicates the normalized magnitude of the rate-of-deformation tensor at the fiber surfaces. At a linear strain of 0.22 [Fig. 10(a)], the fibers are still randomly oriented, and the magnitude of the rate-of-deformation tensor at the fibers' tips is about 3–4 times the magnitude of the rate-of-deformation tensor in case no fibers are present in the matrix. As the strain increases, the fibers orient, thereby increasing the magnitude of the rate-of-deformation tensor at the fibers' tips. For oriented fibers, see Fig. 10(e), the magnitude of the rate-of-deformation tensor obtains its minimum at the center of the fibers due to symmetry of the flow and increases toward the fibers' tips, which is mainly contributed by the *D _{xz}*,

*D*,

_{zx}*D*, and

_{yz}*D*components. The local shear rates increase due to the no-slip boundary condition at the fibers' surfaces. The location of the fibers is also particularly interesting, because at higher strains a microstructural change is visible, where one can see regions without fibers and regions with fiber clusters.

_{yz}Corresponding to Fig. 10, in Fig. 11, the trace of the viscoelastic stress tensor is shown normalized by $ \tau \xaf$, where $ \tau \xaf$ describes the time-dependent trace of the viscoelastic stress tensor in case no fibers are present in the matrix. In this figure, regions under extension have a normalized stress above one, and regions under compression have a normalized stress below one. The snapshots are 2 D slices in the *xy* plane, showing the normalized stress in the matrix. For clarity and interpretation reasons, also the 3 D fibers close to the sliced plane are shown. High extensional stresses occur between the tips of passing fibers, as is visible on the right-hand side of Fig. 11(d). The high extensional stress spreads over a larger region when fibers separate as is visible in Fig. 11(e). In the contraction direction, perpendicular to the stretch direction, fibers approach, which results in compressive stresses in the matrix.

### C. Fiber orientation kinetics

*L*/

*D*) of 10 with random initial positions and orientations. Hencky strain rates between 0.1 and 10 s

^{−1}are considered, and fiber volume fractions are varied between 1% and 10%. The initial domain size in these simulations is fixed, according to

^{1}Moreover, the black solid line in the figure represents the orientation kinetics described by a multi-particle model. In this model, every single fiber follows the Jeffery motion, and fiber–fiber interactions are neglected, and for details, we refer the reader to our previous work.

^{1}

First, the influence of the Weissenberg number on the fiber orientation kinetics is investigated for an LDPE matrix with fiber volume fraction of 1% in Fig. 12(a). Weissenberg numbers between 5.3 and 530 are considered by applying Hencky strain rates of 0.1–10 s^{−1}. The figure shows that all solid lines in Fig. 12(a) overlap, indicating that there is no influence of the Weissenberg number on the fiber orientation kinetics. Moreover, the simulation results show excellent agreement with the experiments.

Second, the influence of strain hardening on the fiber orientation kinetics is investigated in Fig. 12(b). These simulations consider a fiber volume fraction of 1% and a Hencky strain rate of 1 s^{−1} with either an LDPE or an LLDPE matrix. The numerically obtained fiber orientation kinetics (solid lines) in the figure overlap, implying that strain hardening of the matrix does not significantly alter the orientation kinetics. Again the simulation results show excellent agreement with the experiments.

Finally, the influence of fiber volume fraction on the fiber orientation kinetics is investigated using a Hencky strain rate of 1 s^{−1} in Fig. 12(c). The maximum linear strain in the simulations with 10% fiber volume fraction is limited to about 3–4, due to computational difficulties at high volume fractions. For the simulations with a fiber volume fraction of 1%, the maximum linear strain is about 10. Up until a linear strain of 3–4, the fiber orientation kinetics of simulations with 1% or 10% fiber volume fraction overlap, implying that there is no influence of fiber volume fraction on the fiber orientation kinetics for the range studied here.

To conclude this section, no influence of the Weissenberg number, strain hardening, and fiber volume fraction on the fiber orientation kinetics is visible in the simulations, which is in line with our previous experimental study.

### D. Bulk rheology

In this section, the bulk rheology of fiber-filled polymers is investigated. First, direct numerical simulations with a Newtonian matrix are considered and the resulting bulk rheology is compared to the bulk rheology of Hinch and Leal's constitutive equation for fiber composites, see Sec. V. Then, a viscoelastic matrix is considered and the bulk rheology is investigated for the same simulations as performed in Sec. VI C and these numerical results are compared to experimental results obtained in our previous study.^{1}

#### 1. Newtonian matrix

First, the bulk rheology is investigated for fibers in a Newtonian matrix. To investigate the contribution of initial orientation to the bulk rheology, the initial fiber orientation is controlled with random *ψ*_{0} and fixed *θ*_{0} orientation angles. Furthermore, the fiber volume fraction is 10%, the fiber aspect ratio (*L*/*D*) is 10, and the number of fibers in the simulation is 24. Moreover, a Hencky strain rate of 1 s^{−1} is used and the Newtonian viscosity equals 19 500 Pa s, which corresponds to the zero shear viscosity of LLDPE. The transient uniaxial extensional viscosity is shown as a function of *θ*_{0} by the solid lines of Figs. 13(a) and 13(b). The FEM simulations predict the highest viscosity when the fibers are aligned in the stretch direction. Also, when fibers are aligned perpendicular to the stretch direction, they significantly increase the viscosity by their resistance against contractional flow. The minimum viscosity is obtained when the fibers' resistance against flow in the stretch direction and contraction direction is minimal, which is when the fibers are oriented with $ \theta = cos \u2009 ( 1 / 3 ) \u2248 54.7$°. This is clearly visible as a dip in viscosity for fibers that start with $ \theta 0 > 54.7$°. Up to 1 s, the strain hardening in FEM simulations is purely a result of fiber orientation kinetics.

The bulk rheology obtained from the FEM simulations is compared to the bulk rheology resulting from the constitutive equation of Sec. V. The parameters of this constitutive equation ( $ N i , \u2009 N p$, and $ N s$) have been derived by Dinh and Armstrong^{26} for a slender body approach in the semi-concentrated regime. They concluded that for a fiber aspect ratio (*L*/*D*) of 10 and a fiber volume fraction ( $\varphi $) of 10%, $ N p \varphi = 1.93$, and $ N s \varphi \u226a N p \varphi $, which means that the $ N s$ term can be neglected.^{17,54} Furthermore, the slender-body approach assumes the thickness of the fibers to be zero, as a result $ N i = 0$. The resulting transient uniaxial extensional viscosity for different initial fiber orientations is shown by the dotted lines of Fig. 13(a). Here, the initial fiber orientation tensor is computed using Eq. (25) with 1000 fibers pseudorandomly oriented using a random *ψ*_{0} and a fixed *θ*_{0} as shown in the legend of the figure. The FEM simulations and the Dinh and Armstrong model show the same qualitative behavior; however, the slender-body model of Dinh and Armstrong under-predicts the initial viscosity because it does not take the isotropic contribution of the fibers to the viscosity into account ( $ N i = 0$).

The constitutive model of Eq. (55) is improved by fitting $ N i$ and $ N p$ to the FEM simulations with initial fiber orientations of $ \theta 0 = 60$° and $ \theta 0 = 0$°, respectively. The resulting transient uniaxial extensional viscosity for different initial fiber orientations is shown by the dashed line in Fig. 13(b), where $ N i = 2$ and $ N p = 35$. The model agrees very well with the FEM simulations. However, at a time of about 1 s, additional strain hardening is visible in the FEM simulations, even though all fibers are aligned in the stretch direction, whereas the constitutive model predicts a plateau in the viscosity.

To investigate this effect further, especially to isolate finite size effects in the simulations, in Fig. 13(c) the normalized transient uniaxial extensional viscosity is shown for Newtonian FEM simulations with different numbers of fibers in the RVE. All fibers are initially oriented in the stretch direction ( $ \theta 0 = 0$°). For increasing number of fibers in the simulation, the strain hardening occurs at the same time, and with the same slope, which indicates that strain hardening is not caused by finite size effects. At the time of the strain hardening, a microstructural change is visible whereby fibers cluster together, as discussed in Sec. VI B. This local densification, whereby fiber clusters act as a dispersed phase with a larger effective volume fraction, due to entrapped matrix material, increases the resistance against flow and can potentially contribute to the strain hardening. Another explanation might be that in between approaching fibers in the contraction direction, perpendicular to the stretch direction, high local rates occur. These high local rates exponentially increase with smaller distances between neighboring fibers, and result in higher stresses. It is, however, also possible that the amount of fibers in the simulations is still insufficient to describe the rheology correctly, as the triperiodic boundary conditions prescribe that the fiber clusters span the full cross-sectional area at the macroscopic level. At other locations, the triperiodic boundary conditions prescribe no presence of fibers across the cross-sectional area. In a real system, the appearance of many local fiber clusters at different locations is more likely. Increasing the number of fibers in the simulations even further was due to computational limitations not possible.

In Fig. 14, the transient uniaxial extensional viscosity for different initial fiber orientations is investigated for various fiber volume fractions, $\varphi $, and fiber aspect ratios, *L*/*D*. Here, single Newtonian direct numerical FEM simulations (solid lines) agree well with the constitutive model of Eq. (55), with $ N i = 2 , \u2009 N p = ( L / D ) 1.55$, and $ N s = 0$ (dashed lines), implying that up to the concentrated regime, $ \varphi \u2248 D / L$, the Tuckers constitutive model with these parameters can be used to describe the orientation-induced transient uniaxial extensional viscosity of fiber-filled composites.

#### 2. Viscoelastic matrix

In this section, the transient uniaxial extensional viscosity is investigated for the same simulations as performed in Sec. VI C. Figures 15 and 16 show the numerically obtained transient uniaxial extensional viscosity (solid lines) and the experimentally obtained transient uniaxial extensional viscosity (markers) for LDPE composites and LLDPE composites, respectively. The experimental data are obtained in our previous work.^{1} In Fig. 15(a), the transient uniaxial extensional viscosity of LDPE with a fiber volume fraction of 1% is shown for different Hencky strain rates. The experimental data are well captured by the numerical simulations. Furthermore, in the simulation the transient uniaxial extensional viscosity is generated for a Hencky rate of 10 s^{−1}, which cannot be obtained experimentally but is relevant for processing. For simulations performed with a Hencky strain rate of 0.1 s, a time step of 0.1 s is used, which prevents a direct comparison with the experimental data up to approximately 0.5 s. The data in this range are therefore indicated with a dashed line. Figure 15(b) shows the transient uniaxial extensional viscosity of the suspension obtained with a Hencky rate of 1 s^{−1} for LDPE with different fiber volume fractions. Also in this figure, the experimental and numerical data agrees well. Moreover, Fig. 15(b) shows that strain hardening due to fiber orientation is negligible compared to the strain hardening behavior of the unfilled viscoelastic matrix.

In Fig. 16, the transient uniaxial extensional viscosity is shown for the non-strain-hardening matrix LLDPE with different fiber volume fractions. In these results, the applied Hencky rate is 1 s^{−1} with a random initial fiber orientation. The figure shows that the numerical viscoelastic simulations under-predict the experimentally obtained transient uniaxial extensional viscosity. This originates from under-prediction of the Giesekus material model on the pure LLDPE. The simulations with a fiber volume fraction of 10% show a significant amount of strain hardening, which is mostly related to changes in fiber orientation: from a random fiber orientation to an oriented state in the stretch direction. At a time of about 1 s, the strain hardening can also be attributed to the local change in microstructure.

## VII. CONCLUSION AND DISCUSSION

In this study, a numerical framework was developed and used to investigate the fiber orientation kinetics and the resulting rheology of fiber-filled polymer composites in uniaxial extension. The fiber orientation kinetics are not influenced by the Weissenberg number ( $ W i = 0.53 \u2212 53$), strain hardening, and fiber volume fraction ( $ \varphi = 1 % \u2212 10 %$), which agrees very well with experimental results obtained in our previous study.^{1} As a consequence, in uniaxial extension and for the concentration regime considered in this study, the phenomenological terms slowing down the orientation kinetics in the Folgar–Tucker models are not required, which is in line with the observations of Lambert and Baird,^{55} who investigated orientation kinetics in planar extension far in the concentrated regime ( $ \varphi \u226b D / L$). High agreement was also obtained between the rheology predicted by the numerical framework and the experimental results obtained in our previous study,^{1} except for systems with a high fiber volume fraction and a non-strain-hardening matrix. Here, in the experiments, stress localization occurred, resulting in local sample failure at earlier strains.^{5} In contrast, the numerical framework assumed perfect adhesion between the fibers and the matrix, and predicted orientation-induced strain hardening. Orientation-induced strain hardening was, up to the concentrated regime ( $ \varphi \u2248 D / L$), accurately described by Hinch and Leal's constitutive model with parameters predicted by the numerical framework. One crucial predicted parameter is the isotropic contribution of the fibers to the viscosity. At large strain for aligned fibers, the numerical framework predicts additional strain hardening resulting from a local change in microstructure as fiber clusters start appearing, e.g., local densification. This effect might be more pronounced in highly filled polymer fiber composites and is not captured in the current rheological fiber models.

The numerical framework is also able to determine the time-dependent conformation tensor per polymer mode, which, for semi-crystalline matrices, could give insights into crystallization kinetics. As a result, the numerical framework presented in this study could be used to investigate the complex interplay between flow, fiber orientation, rheology, and crystallization kinetics, and could thereby yield an accurate prediction of the microstructure of fiber-filled polymer composites.

## SUPPLEMENTARY MATERIAL

See the supplementary material for additional numerical results, comparing different closure approximations as well as a multi-particle model.

## ACKNOWLEDGMENTS

This research forms part of the research program of DPI, Project No. 840 Engineering the rheology ANd processinG-induced structural anisotropy of poLymEr composites with non-Brownian fibrous particles (ANGLE). The authors thank M. A. Hulsen at the Eindhoven University of Technology (TU/e) for access to the TFEM software libraries. Furthermore, this work used the Dutch national e-infrastructure with the support of the SURF Cooperative using Grant No. EINF-5641.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Thijs R. N. Egelmeers:** Conceptualization (supporting); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (supporting). **Ruth Cardinaels:** Conceptualization (equal); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). **Patrick D. Anderson:** Conceptualization (supporting); Formal analysis (supporting); Funding acquisition (supporting); Investigation (supporting); Methodology (supporting); Project administration (equal); Resources (equal); Software (supporting); Supervision (supporting); Writing – review & editing (supporting). **Nick O. Jaensson:** Conceptualization (equal); Formal analysis (equal); Funding acquisition (supporting); Investigation (equal); Methodology (equal); Project administration (equal); Resources (equal); Software (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: DETERMINATION OF THE RELATIVE RVE SIZE

In this section, viscoelastic uniaxial extensional simulations are performed to determine the influence of RVE size, expressed as the number of particles in the domain, on the averages and standard deviations of the fiber orientation and the bulk stresses. In these simulations, spherocylindrical particles with aspect ratio, *L*/*D*, of 10 are used with a fiber volume fraction of 5%. Moreover, a single-mode Giesekus material model with $ \eta p = 10 5$ Pa s, *λ* = 1 s, $ \alpha = 0.005 , \u2009 \eta sol = 0$ Pa s is used, and a Hencky strain rate of $ \epsilon \u0307 H = 1$ s^{−1} is applied. RVEs with 8, 16, and 32 particles were investigated, which fixes the initial RVE size, $ L 0 , x$, according to Eq. (59). Eight particles are the lowest amount of particles for a fiber volume fraction of 5% since lowering the amount of particles resulted in domain sizes smaller than the length of a single particle. The positions and orientations of the particles are randomly generated. The average (solid line) and standard deviation on the average (dashed line) for the anisotropy are shown in Fig. 17(a) and the average (solid line) and standard deviation on the average (dashed line, not visible due to the log scale) for the transient uniaxial extensional viscosity is shown in Fig. 17(b), which are computed using 12 realizations. The standard deviation of the anisotropy decreases with applied strain. At a strain of 2, these are equal to 0.032, 0.0176, and 0.014 for RVE's with, respectively, 8, 16, and 32 particles. The standard deviation can be decreased by either increasing the number of realizations or by increasing the amount of particles in the RVE.

## REFERENCES

*Fundamentals of Fiber Orientation: Description, Measurement and Prediction*

*Rheology: Principles, Measurements, and Applications*

*The Art of Molecular Dynamics Simulation*