Numerical simulation of fiber orientation kinetics and rheology of fiber-filled polymers in uniaxial extension

Egelmeers, Thijs R. N.; Cardinaels, Ruth; Anderson, Patrick D.; Jaensson, Nick O.

doi:10.1063/5.0189644

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,¹⁶ 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.¹⁷ 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.¹⁸

As Kruger et al.¹⁸ 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¹⁹ 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⁵ 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,²⁴ 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²⁵ 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.,¹ Wagner et al.,³² and Wang et al.³³ did not show the orientation-induced strain hardening. In contrast with these findings, Ferec et al.³⁰ 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⁵ 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.¹ The rheology resulting from the FEM simulations will be compared to the rheology from a constitutive model proposed by Hinch and Leal.² 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

A suspension of rigid non-Brownian spherocylindrical or prolate spheroidal particles is considered with a solid volume fraction

ϕ

⁠. The suspension is subjected to a macroscopic uniaxial extensional flow field, as shown in Fig. 1. The flow will introduce stress in the matrix, and particles will rotate and translate accordingly. In this work, a triperiodic domain with size

L_{x} (t), L_{y} (t)

⁠, and

L_{z} (t)

is considered at the micro-scale assuming that the macroscopic length scale is much larger than the microscopic length scale. The velocity in the considered domain (u) is equal to

u^{*} + \hat{u}

⁠, where

u^{*}

is the periodic part of the velocity, and

\hat{u}

the velocities described by macroscopic uniaxial extensional field, according to

\hat{u} (x, y, z) = {[- x \frac{{\dot{ε}}_{H}}{2}, - y \frac{{\dot{ε}}_{H}}{2}, z {\dot{ε}}_{H}]}^{T},

(1)

where

{\dot{ε}}_{H}

is the Hencky strain rate, and 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

Δ u (x, y, z) = {[- i \frac{{\dot{ε}}_{H}}{2} L_{x} (t), - j \frac{{\dot{ε}}_{H}}{2} L_{y} (t), k {\dot{ε}}_{H} L_{z} (t)]}^{T},

(2)

where

L_{x} (t), L_{y} (t)

⁠, and

L_{z} (t)

are the triperiodic domain sizes, which are a function of time to correctly account for convection, according to^9,34

\begin{matrix} L_{x} (t) = L_{x, 0} exp (\frac{- {\dot{ε}}_{H}}{2} t), \\ L_{y} (t) = L_{y, 0} exp (\frac{- {\dot{ε}}_{H}}{2} t), \\ L_{z} (t) = L_{z, 0} exp ({\dot{ε}}_{H} t), \end{matrix}

(3)

where

L_{x, 0}, L_{y, 0}

⁠, and

L_{z, 0}

are the initial domain sizes. In the remainder of the paper, L_x, L_y, and L_z will be used to indicate the time-dependent domain sizes.

FIG. 1.

View large Download slide

The macroscopic uniaxial extensional field (left) and the microscopic triperiodic domain (right).

A. Particles in the triperiodic domain

A triperiodic domain with particles is considered. The orientation of a particle is described using

p_{i} = {[sin (ψ) sin (θ), cos (ψ) sin (θ), cos (θ)]}^{T},

(4)

where

p_{i}

is the orientation vector for the ith 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. .

FIG. 2.

View large Download slide

Schematic representation of the triperiodic domain. (a) Important definitions used to describe the geometry, the fluid domain is denoted with Ω, and the volume occupied by the ith particle with F_i, (b) continuity of the triperiodic domain with particles separated by the boundary of the triperiodic domain.

B. Balance equations

The mass and momentum balance are used to describe the flow of the fluid. For an incompressible fluid, neglecting inertia and assuming no body forces, the balance equations yield

\nabla \cdot u = 0 in Ω,

(5)

\nabla \cdot σ = 0 in Ω,

(6)

where u is the fluid velocity vector and

σ

is the Cauchy stress tensor

σ = - p I + 2 η_{sol} D + τ,

(7)

where p is the pressure, I the unit tensor,

η_{sol}

the Newtonian solvent viscosity, and

D = (\nabla u + {(\nabla u)}^{T}) / 2

represents the rate of deformation tensor. The term

2 η_{sol} D

is the Newtonian stress tensor and

τ

is the viscoelastic stress tensor. The viscoelastic stress tensor for the Giesekus constitutive model in conformation tensor representation is given by³⁵

τ = \frac{η_{p}}{λ} (c - I),

(8)

λ \overset{▽}{c} + c - I + α {(c - I)}^{2} = 0,

(9)

where c is the conformation tensor,

η_{p}

the polymer viscosity, and λ is the polymer relaxation time. The latter two are related through the polymer modulus

G = η_{p} / λ

⁠. Furthermore, α is a (dimensionless) constitutive parameter, which influences shear thinning behavior and the symbol (⁠

^{▽}

⁠) denotes the upper-convected time derivative

\overset{▽}{c} = \frac{\partial c}{\partial t} + u \cdot \nabla c - {(\nabla u)}^{T} \cdot c - c \cdot (\nabla u),

(10)

where

\partial () / \partial t

is the partial time derivative at a fixed location in space. In the remainder of the paper, the conformation tensor evolution equation [Eq. ] is written as

\frac{\partial c}{\partial t} + u \cdot \nabla c - g (G, c) = 0

with

G = {(\nabla u)}^{T}

⁠. Finally, the pressure level is set using the following constraint:

\int_{Ω} p d V = 0.

(11)

C. Boundary and initial conditions

The balance equations are closed using the boundary conditions and initial condition described in this section. First, triperiodic boundary conditions are imposed on the external domain, which can be expressed as (see Fig. 2)

\begin{matrix} u_{x} |_{Γ_{4}} = u_{x} |_{Γ_{2}} + {\dot{ε}}_{H} L_{x} / 2, \\ u_{y} |_{Γ_{4}} = u_{y} |_{Γ_{2}}, \\ u_{z} |_{Γ_{4}} = u_{z} |_{Γ_{2}}, \end{matrix}

(12)

\begin{matrix} u_{x} |_{Γ_{1}} = u_{x} |_{Γ_{3}}, \\ u_{y} |_{Γ_{1}} = u_{y} |_{Γ_{3}} + {\dot{ε}}_{H} L_{y} / 2, \\ u_{z} |_{Γ_{1}} = u_{z} |_{Γ_{3}}, \end{matrix}

(13)

\begin{matrix} u_{x} |_{Γ_{6}} = u_{x} |_{Γ_{5}}, \\ u_{y} |_{Γ_{6}} = u_{y} |_{Γ_{5}}, \\ u_{z} |_{Γ_{6}} = u_{z} |_{Γ_{5}} - {\dot{ε}}_{H} L_{z}, \end{matrix}

(14)

for the velocities. For the traction vector

t = σ \cdot n

⁠, with n the outward pointed normal on surface

Γ_{i}

⁠, continuity of the stress yields

\begin{matrix} t |_{Γ_{4}} = - t |_{Γ_{2}}, \\ t |_{Γ_{1}} = - t |_{Γ_{3}}, \\ t |_{Γ_{6}} = - t |_{Γ_{5}}, \end{matrix}

(15)

and for the conformation tensor, we obtain

\begin{matrix} c |_{Γ_{4}} = c |_{Γ_{2}}, \\ c |_{Γ_{1}} = c |_{Γ_{3}}, \\ c |_{Γ_{6}} = c |_{Γ_{5}} . \end{matrix}

(16)

The absolute level of the velocities still needs to be set. One way is by prescribing the velocity in a single point in the domain, as was performed in Ref. 36. Another approach is setting the translational velocity of one particle to zero, as described by Jaensson et al.⁹ In this work, the latter approach has been chosen.

It is assumed that the particle motion is slow so that inertia can be neglected. The resulting total hydrodynamic force and torque acting on the complete particle will therefore be zero, respectively,

\int_{\partial F_{i}} σ \cdot n_{i} d S = 0 on \partial F_{i},

(17)

\int_{\partial F_{i}} (x - X_{i}) \times (σ \cdot n_{i}) d S = 0 on \partial F_{i},

(18)

where

\partial F_{i}

is the ith particle surface,

n_{i}

the outwardly directed unit normal vector on that particle surface, x are coordinates on that particle surface, and

X_{i} = {[X_{i}, Y_{i}, Z_{i}]}^{T}

is the local center point coordinate of particle i. Assuming no-slip, the boundary condition for the velocity on the complete particle boundary

\partial F_{i}

will therefore become

u = U_{i} + ω_{i} \times (x - X_{i}) on \partial F_{i},

(19)

where

ω_{i}

is the unknown angular velocity of the ith 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. and . Moreover, in this case, the velocity difference for different periodic center points,

U_{i}

⁠, is corrected using Eq. .

Finally, a stress-free initial condition is assumed, i.e., at the initial time, the viscoelastic stress is zero everywhere in the fluid

τ |_{t = 0} = 0 in Ω .

(20)

D. Arbitrary Lagrange–Euler (ALE) formulation

The domain is described with a boundary-fitted mesh that is moved in time in such a way that the mesh moves with the particles, but not necessarily with the fluid. Therefore, the governing equations have to be rewritten in the Arbitrary Lagrange–Euler (ALE) formulation.³⁷ 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

\overset{▽}{c} = {\frac{\partial c}{\partial t} |}_{x_{m}} + (u - u_{m}) \cdot \nabla c - {(\nabla u)}^{T} \cdot c - c \cdot (\nabla u),

(21)

where the first term on the right-hand side is the partial derivative for a constant mesh coordinate

x_{m}

⁠, and

u_{m}

is the mesh velocity.

E. Bulk stresses

The bulk stress is the volume-averaged stress in the domain and is computed using Batchelor's formula, according to³⁸

⟨ σ ⟩ = \frac{1}{V} \int_{Ω \cup F} σ d V = \frac{1}{V} \int_{Ω} σ d V + \frac{1}{V} \int_{F} σ d V,

(22)

where the brackets

⟨ . ⟩

indicate a volume-averaged quantity, V is the volume of the domain, Ω is the volume occupied by the fluid in the domain, and

F = F_{1} \cup F_{2} \cup \dots F_{N}

is the volume occupied by the particles in the domain, see Fig. 2(a). By substituting

σ = \nabla \cdot (σ^{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

\partial F

⁠. Furthermore, Eq. is substituted in Eq. ,^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

⟨ σ ⟩ = - p' I + 2 η_{sol} D' + τ' + \sum_{i = 1}^{N} \int_{\partial F_{i}} σ \cdot n (x - X_{i}) d S,

(23)

where the prime denotes the integral over the fluid volume divided by the entire volume, e.g.,

D' = \int_{Ω} D d V / V

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

⟨ η_{E} ⟩ = \frac{⟨ σ_{z z} ⟩ - ⟨ σ_{x x} ⟩}{{\dot{ε}}_{H}},

(24)

where

⟨ η_{E} ⟩

is the uniaxial extensional viscosity.

F. Anisotropy

The anisotropy is defined as the average level of particle alignment in the stretching direction, analogous to Egelmeers et al.¹ The average particle orientation is computed using the particle orientation tensor¹⁰

P = \frac{1}{N} \sum_{i = 1}^{N} p_{i} p_{i} .

(25)

For extension along the z axis with initial random particle orientation, all off-diagonal components of P 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 x and y directions with the average orientation component in the z direction

ξ_{s} = 1 - \frac{\sqrt{ν_{x}^{2} + ν_{y}^{2}}}{ν_{z}} \approx 1 - \frac{\sqrt{P_{x x}^{2} + P_{y y}^{2}}}{P_{z z}},

(26)

where

ν_{x}, ν_{y}

⁠, and ν_z 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,

ξ_{s} = - 0.41

and when all particles are perfectly aligned in the stretch direction,

P_{x x} = P_{y y} = 0

and P_zz = 1, so that

ξ_{s} = 1

⁠.

G. RVE statistics

In this section, the criteria for the number of realizations for each RVE (Representative Volume Element) per parameter set are defined to obtain meaningful average quantities for the particle orientation and the bulk stress evolution. The standard deviation on the average for the anisotropy and the bulk stress is defined as⁹

s_{{\bar{ξ}}_{s}} (t) = \frac{1}{\sqrt{n_{r}}} \sqrt{\frac{1}{n_{r} - 1} \sum_{i = 1}^{n_{r}} {(ξ_{s} (t) |_{i} - {\bar{ξ}}_{s} (t))}^{2}},

(27)

s_{{\bar{η}}_{E}^{+}} (t) = \frac{1}{\sqrt{n_{r}}} \sqrt{\frac{1}{n_{r} - 1} \sum_{i = 1}^{n_{r}} {(η_{E}^{+} (t) |_{i} - {\bar{η}}_{E}^{+} (t))}^{2}},

(28)

respectively, where

n_{r}

is the number of realizations,

{\bar{ξ}}_{s}

and

{\bar{η}}_{E}^{+}

are the average anisotropy and the average transient uniaxial extensional viscosity, and

ξ_{s} (t) |_{i}

and

η_{E}^{+} (t) |_{i}

are the anisotropy and transient uniaxial extensional viscosity for the ith realization. The interval that contains the true average with 95% confidence equals

{\bar{ξ}}_{s} \pm 2 s_{{\bar{ξ}}_{s}}

and

{\bar{η}}_{E}^{+} \pm 2 s_{{\bar{η}}_{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:

e_{ξ_{s}} = 2 s_{{\bar{ξ}}_{s}} < 0.04,

(29)

e_{η_{E}^{+}} = \frac{2 s_{{\bar{η}}_{E}^{+}}}{{\bar{η}}_{E}^{+}} < 0.04,

(30)

where

e_{ξ_{s}}

and

e_{η_{E}^{+}}

are the relative errors of the anisotropy and the transient uniaxial extensional viscosity, respectively. In the remainder of the paper, the overbar notation will be omitted to denote average values and only the averages will be considered, unless otherwise indicated.

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

The weak formulation is obtained by multiplying the balance equations with a test function and integrating over the domain. After using Gauss' theorem and integrating by parts, the weak formulation can be expressed as follows: find

p, u, U_{i}, ω_{i}, c

⁠, and

λ_{i}

(where

i = 1, \dots, N

⁠), such that

(q, \nabla \cdot u) = 0,

(31)

\begin{matrix} ({(\nabla v)}^{T}, ν (\nabla u - G^{T})) + (D_{v}, 2 η_{sol} D + τ) \\ - (\nabla \cdot v, p) + {⟨ v - (χ_{1} \times (x - X_{1})), λ_{1} ⟩}_{\partial F_{1}} \\ + {⟨ v - (V_{2. . N} + χ_{2. . N} \times (x - X_{2. . N})), λ_{2. . N} ⟩}_{\partial F_{2. . N}} = 0, \end{matrix}

(32)

(H, - \nabla u + G^{T}) = 0,

(33)

\begin{matrix} (d + τ (u - u_{m}) \cdot \nabla d, {\frac{\partial s}{\partial t} |}_{x_{m}} + (u - u_{m}) \cdot \nabla s - h (G, s)) = 0, \end{matrix}

(34)

{⟨ μ_{1}, u - ω_{1} \times (x - X_{1}) ⟩}_{\partial F_{1}} = 0,

(35)

{⟨ μ_{2. . N}, u - (U_{2. . N} + ω_{2. . N} \times (x - X_{2. . N})) ⟩}_{\partial F_{2. . N}} = 0,

(36)

for all test functions, q,

v

⁠,

V_{i}, μ_{i}

d,

χ_{i}

⁠, H. Moreover,

G = \nabla u^{T}

is the velocity gradient,

D_{v} = \frac{1}{2} (\nabla v + {(\nabla v)}^{T}), (\cdot, \cdot)

denotes the inner product in domain Ω,

{⟨ \cdot, \cdot ⟩}_{\partial F}

denotes the inner product on particle surface

\partial F

⁠, τ and ν are parameters for the SUPG and DEVSS-G stabilization,

s = log c

⁠, and

h (G, s)

is a tensor function.⁴³

B. Spatial discretization

Isoparametric tetrahedral Taylor–Hood elements are used for the velocity and the pressure, employing quadratic interpolation of the velocities (P₂) and linear interpolation for the pressures (P₁). For the velocity gradient and the conformation tensor, P₁ interpolation is used. The mesh is generated using Gmsh.⁴⁵ 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

In this section, the time stepping procedure will be discussed. An implicit stress formulation and a second-order semi-implicit Gear scheme are used for the conformation prediction, as described by D'Avino and Hulsen.⁴⁹ In addition, to not be limited in time step size by very fast modes, the terms

s^{n + 1} / λ

and

{\hat{s}}^{n + 1} / λ

are added on the left and right sides of the equation, respectively. The second-order semi-implicit Gear scheme in log form now becomes

\begin{matrix} \frac{3}{2} \frac{s^{n + 1}}{Δ t} + u^{n + 1} \cdot \nabla s^{n + 1} + \frac{s^{n + 1}}{λ} = \frac{2 s^{n} - (1 / 2) s^{n - 1}}{Δ t} + h (G^{n + 1}, {\hat{s}}^{n + 1}) + \frac{{\hat{s}}^{n + 1}}{λ}, \end{matrix}

(37)

where the current time step is indicated with n and the value of a variable at the current time step with a superscript notation: i.e.,

u |_{t = n Δ t} = u^{n}

⁠. Moreover,

{\hat{s}}^{n + 1} = 2 s^{n} - s^{n - 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,

λ | D |

⁠.

Quaternions are used to describe the orientation of the particles. Quaternions are an extension of complex numbers, which have one real axis and three imaginary axes, which can be written as

q = {[q_{r}, q_{1}, q_{2}, q_{3}]}^{T}

⁠. The relation between the time derivative of a quaternion and the angular velocity of a particle is given by⁵⁰

{\dot{q}}_{i} = \frac{\partial q_{i}}{\partial t} = \frac{1}{2} W (q_{i}) ω_{i},

(38)

with W a matrix given by

W (q) = [\begin{matrix} - q_{1} & - q_{2} & - q_{3} \\ q_{r} & q_{3} & - q_{2} \\ - q_{3} & q_{r} & q_{1} \\ q_{2} & - q_{1} & q_{r} \end{matrix}],

(39)

where subscript 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

\frac{q_{i}^{1}}{Δ t} = \frac{q_{i}^{0}}{Δ t} + \frac{1}{2} [W (q_{i}^{0}) ω_{i}^{0}] .

(40)

From the second time step onward, a second-order Adams–Bashforth time integration scheme is used, according to

\frac{q_{i}^{n + 1}}{Δ t} = \frac{q_{i}^{n}}{Δ t} + \frac{1}{2} [\frac{3}{2} W (q_{i}^{n}) ω_{i}^{n} - \frac{1}{2} W (q_{i}^{n - 1}) ω_{i}^{n - 1}] .

(41)

After normalizing the quaternion, the rotation matrix,

R_{i}^{n + 1 \leftarrow n}

⁠, is constructed, which maps the particle orientation at the current time step, n, to the particle orientation at the new time step, n + 1, for the ith 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^{- 1} = R^{T}

⁠), and subsequently, the initial state is mapped to the new time step, according to

R_{i}^{n + 1 \leftarrow n} = R_{i}^{n + 1 \leftarrow 0} \cdot R_{i}^{0 \leftarrow n} = R (q_{i}^{n + 1}) \cdot R^{T} (q_{i}^{n}),

(42)

with⁵¹

R (q) = 2 \cdot [\begin{matrix} \frac{1}{2} - q_{2}^{2} - q_{3}^{2} & q_{1} q_{2} - q_{r} q_{3} & q_{1} q_{3} + q_{r} q_{2} \\ q_{1} q_{2} + q_{r} q_{3} & \frac{1}{2} - q_{1}^{2} - q_{3}^{2} & q_{2} q_{3} - q_{r} q_{1} \\ q_{1} q_{3} - q_{r} q_{2} & q_{2} q_{3} + q_{r} q_{1} & \frac{1}{2} - q_{1}^{2} - q_{2}^{2} \end{matrix}],

(43)

where

R (q)

is the rotation matrix that describes the mapping of the initial particle orientation to the particle orientation described by q. The orientation of the particle is updated according to

p_{i}^{n + 1} = R_{i}^{n + 1 \leftarrow n} \cdot p_{i}^{n} .

(44)

The position of the particles' center is updated using a second-order Adams–Bashforth time integration scheme

\frac{X_{i}^{n + 1}}{Δ t} = \frac{X_{i}^{n}}{Δ t} + (\frac{3}{2} U_{i}^{n} - \frac{1}{2} U_{i}^{n - 1}) .

(45)

Subsequently, the coordinates of the nodes on the particles' surface are updated according to

x^{n + 1} = X_{i}^{n + 1} + R_{i}^{n + 1 \leftarrow n} \cdot (x^{n} - X_{i}^{n}) on \partial F_{i} .

(46)

The coordinates of the nodes on the surface of the domain are updated, while keeping the location of particle 1 constant, according to

x^{n + 1} = (x^{n} - X_{1}) \cdot A + X_{1} on Γ_{i},

(47)

where A is a diagonal tensor with

A_{x x} = A_{y y} = exp (- {\dot{ε}}_{H} Δ t / 2)

and

A_{z z} = exp ({\dot{ε}}_{H} Δ t)

⁠. Note that nodes shared by particles and periodic surfaces are updated according to the particle displacement.

Finally, a Laplace problem is solved³⁷

\nabla \cdot (k_{e} \nabla (Δ x_{m})) = 0 in Ω,

(48)

using the nodal displacements at the surfaces of Eqs. and as boundary conditions (⁠

Δ x_{m} = x_{m}^{n + 1} - x_{m}^{n}

⁠). Moreover, k_e 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.

D. Remeshing and projection

As the particles move and rotate, the mesh deforms. Eventually, the elements in the mesh can get too distorted, which results in numerical instabilities. To prevent numerical instabilities from occurring, the mesh distortion is determined every time step by tracking the change in volume and the change in aspect ratio of the individual mesh elements³⁷

f_{1}^{e} = | log (V^{e} / V_{0}^{e}) |,

(49)

f_{2}^{e} = | log (S^{e} / S_{0}^{e}) |,

(50)

where

V^{e}

is the element volume;

V_{0}^{e}

is the undistorted element volume;

S^{e}

is the element aspect ratio, which is defined as

S^{e} = {(L_{max}^{e})}^{3} / V_{0}^{e}

⁠, with

L_{max}^{e}

being the maximum distance between vertices in an element; and

S_{0}^{e}

is the undistorted element aspect ratio. If an element is too distorted (⁠

f_{1}^{e} > 0.1

or

f_{2}^{e} > 0.1

⁠), the mesh is discarded and a remeshing procedure is performed.

Meshing is, as explained in Sec. III B, performed using Gmsh. After meshing, a projection problem⁵² is solved to determine the solution of the new mesh in the previous time step (⁠ $c^{n}, c^{n - 1}, u^{n}, u^{n - 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.

FIG. 3.

View large Download slide

Schematic representation of the remeshing procedure in the case of translated domain vertices. The dotted box represents the domain vertices of the old mesh. Once the elements are too distorted the mesh is discarded, and a new mesh is generated. Then, a projection problem is solved to determine the coordinates and solutions of the new mesh at the previous time step.

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

In this section, mesh and time convergence studies are performed to ensure accurate simulations with cut particles. A prolate spheroidal particle with aspect ratio

L / D = 4

is considered in a Newtonian matrix with a viscosity of 1 Pa s. Moreover, a Hencky strain rate of 1 s⁻¹ is applied in uniaxial extensional flow. The initial orientation vector of the particle equals

p = {[\frac{1}{2} \sqrt{2}, 0, \frac{1}{2} \sqrt{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,

ω_{y, a} = - \frac{3}{4} {\dot{ε}}_{H} \frac{r_{e}^{2} - 1}{r_{e}^{2} + 1} sin (2 θ),

(51)

θ_{a} (t) = arctan [exp (- \frac{3}{2} {\dot{ε}}_{H} \frac{r_{e}^{2} - 1}{r_{e}^{2} + 1} \cdot t) tan (θ_{0})],

(52)

where

r_{e}

is the particle aspect ratio (L/D). The alignment angle θ describes the angle between the stretch direction and the fiber orientation vector (⁠

θ = arccos (p_{z})

⁠). The evolution of the other angle to fully describe the orientation of the particle, ψ, is not of interest since

ψ (t) = 0

⁠.

FIG. 4.

View large Download slide

Mesh and time convergence. (a) Initial position of the considered triperiodic ellipsoidal particle (L = 4) used in the convergence tests with dimensionless domain size (⁠ $L_{0, x} / L$ ⁠) of 2, (b) mesh convergence error as a function of dimensionless element size, and (c) time convergence error as a function of dimensionless time step.

First, mesh size convergence is studied. For this, the 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

‖ e_{M} ‖ = \frac{| ω_{y, FEM} - ω_{y, a} |}{| ω_{y, a} |} .

(53)

Figure 4(b) shows this error as a function of element size,

h = \sqrt{A_{ps} / n_{e, ps}}

⁠, where

A_{ps}

is the surface area of the spheroid and

n_{e, ps}

is the number of elements on the surface of the spheroid. The different lines in the figure correspond to different initial dimensionless domain sizes

L_{0, x} / L

⁠, where 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.

Time convergence is studied for different element sizes at an initial dimensionless domain size,

L_{0, x} / L

⁠, of 20. The time convergence error

‖ e_{T} ‖

is defined as

‖ e_{T} ‖ = \frac{| θ_{FEM} |_{{\dot{ε}}_{H} t = 1} - θ_{a} |_{{\dot{ε}}_{H} t = 1} |}{| θ_{a} |_{{\dot{ε}}_{H} t = 1} |} .

(54)

The error is thus evaluated at a Hencky strain of 1 (linear strain of 1.72). The results are shown in Fig. 4(c). Here,

n_{d}

is a mesh size parameter. It describes the number of nodes on the diameter of the prolate spheroid. A slope of 2 is visible in Fig. 4(c), which is expected for the second order time integration scheme used. A deviation from the slope of 2 is visible when the mesh size error starts to dominate and occurs at lower

{\dot{ε}}_{H} Δ t

for increasing

n_{d}

⁠. A Hencky strain step,

\dot{ε_{H}} Δ t

⁠, of 0.01 is used in the remainder of this paper, unless otherwise indicated.

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 $η_{p} = 10^{5}$ Pa s, λ = 1 s, $α = 0.005$ ⁠, and $η_{sol} = 0$ Pa s, is used in these simulations. Furthermore, a Hencky strain rate of 1 s⁻¹ 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 - 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.

TABLE I.

Mesh size characteristics for the different considered meshes, characterized by the number of nodes on the diameter of the particle (⁠ $n_{d}$ ⁠), and the number of nodes between particles (⁠ $n_{f - f}$ ⁠), leading to different total numbers of elements in the mesh (⁠ $n_{elem}$ ⁠).

Mesh ID	$n_{d}$	$n_{f - 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 - 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 - f}$ ⁠, is clearly visible by comparing the area between the two particles at the right of Figs. 5(a) and 5(e).

FIG. 5.

View large Download slide

Typical meshes corresponding to mesh size characteristic (a) M1, (b) M2, (c) M3, (d) M4, and (e) M5 of Table I, respectively.

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.

FIG. 6.

View large Download slide

Bulk properties for the mesh size characteristics indicated in Table I. (a) Anisotropy of the fibers in the bulk as a function of linear strain, and (b) transient uniaxial extensional viscosity as a function of time.

FIG. 7.

View large Download slide

Trace of the conformation tensor computed on a cross section in the RVE at a linear strain of 4, for different mesh size characteristics of Table I. (a) M1, (b) M2, (c) M3, (d) M4, and (e) M5.

V. CONSTITUTIVE MODEL FOR RHEOLOGY OF FIBER SUSPENSIONS

The bulk stress resulting from the direct numerical simulations as described in Sec. II E will be compared to the fiber-induced stress predicted using a constitutive model. In the constitutive model, first developed by Hinch and Leal in 1973, a Newtonian matrix is assumed and the stress is a function of the fibers' aspect ratio, volume fraction, and orientation, according to²^,

\begin{matrix} τ = 2 η_{sol} [(1 + N_{i} ϕ) D + N_{p} ϕ P : D \\ + N_{s} ϕ (D \cdot P + P \cdot D)], \end{matrix}

(55)

where

ϕ

is the fiber volume fraction,

N_{i}

is the isotropic contribution of the fibers to the viscosity,

N_{p}

is the particle number which describes the resistance against stretching of a fiber, and

N_{s}

is the shear number, which includes the difference in viscosity for fibers aligned in the vorticity direction and fibers aligned in the velocity direction, hence

N_{s} = 0

for extensional flow. Furthermore,

P

is the fourth-order orientation tensor, which must be computed using a closure approximation.¹⁷ 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 P is a function of time. Its evolution is computed using Jeffery's equation in tensor form¹⁷

\dot{P} = W \cdot P - P \cdot W + β (D \cdot P + P \cdot D - 2 D : P),

(56)

where

W = (\nabla u - {(\nabla u)}^{T}) / 2

is the vorticity tensor and β is the particle shape factor, which is related to geometry and aspect ratio. For the particles in this work, β = 1. By combining Eqs. and , 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.¹ 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.¹ 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.¹ 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.

TABLE II.

Giesekus model parameters describing transient experiments for LDPE at 150 °C.

Mode	λ (s)	$η_{p}$ (Pa s)	α
1	$1.00 \times 10^{- 2}$	$9.29 \times 10^{2}$	0.49
2	$8.51 \times 10^{- 2}$	$3.11 \times 10^{3}$	0.49
3	$5.18 \times 10^{- 1}$	$8.91 \times 10^{3}$	0.49
4	$3.02$	$2.18 \times 10^{4}$	0.25
5	$1.65 \times 10^{1}$	$3.82 \times 10^{4}$	0.11
6	$7.53 \times 10^{1}$	$3.76 \times 10^{4}$	0.11
7	$2.74 \times 10^{2}$	$2.34 \times 10^{4}$	0.08

TABLE III.

Giesekus model parameters describing transient experiments for LLDPE at 150 °C.

Mode	λ (s)	$η_{p}$ (Pa s)	α
1	$1.15 \times 10^{- 2}$	$3.97 \times 10^{3}$	0.49
2	$6.35 \times 10^{- 2}$	$7.11 \times 10^{3}$	0.49
3	$3.18 \times 10^{- 1}$	$5.34 \times 10^{3}$	0.49
4	$2.06$	$1.61 \times 10^{3}$	0.49
5	$1.70 \times 10^{1}$	$5.44 \times 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⁻¹, 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.¹ 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.

FIG. 8.

View large Download slide

Transient extensional viscosity (a), transient shear viscosity (b), and transient first normal stress difference (c) for unfilled LDPE (markers) compared to the corresponding 7-mode Giesekus model (solid lines), (d) shows the complex viscosity as a function of frequency, obtained using small-amplitude oscillatory shear experiments in our previous work¹ (markers) and the steady-state viscosity prediction of the 7-mode Giesekus model, describing LDPE (solid line). Experimental data (markers) are obtained at a temperature of 150 °C.

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.

FIG. 9.

View large Download slide

Transient extensional viscosity (a), transient shear viscosity (b), and transient first normal stress difference (c) for unfilled LLDPE (markers) compared to the corresponding 5-mode Giesekus model (solid line), (d) shows the complex viscosity as a function of frequency, obtained using small-amplitude oscillatory shear experiments in our previous work¹ (markers) and the steady-state viscosity prediction of the 5-mode Giesekus model, describing LLDPE (solid line). Experimental data (markers) are obtained at a temperature of 150 °C.

The characteristic relaxation times and polymer viscosities of LDPE and LLDPE are used to compute the Weissenberg number. The Weissenberg number defines the ratio between elastic and viscous forces in the matrix, according to

Wi = λ_{c} {\dot{ε}}_{H},

(57)

where

λ_{c}

is a characteristic timescale of the matrix. In this study, the viscosity-averaged relaxation time will be used as representative timescale

λ_{c} = \frac{\sum_{i = 1}^{N_{modes}} η_{p, i} λ_{i}}{\sum_{i = 1}^{N_{modes}} η_{p, i}} .

(58)

Using the values of Tables II and III, the viscosity-averaged relaxation times of LDPE and LLDPE are, respectively, 53 and 0.62 s at a temperature of 150 °C.

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⁻¹, 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_yz, and D_yz 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.

FIG. 10.

View large Download slide

Normalized magnitude of the rate-of-deformation tensor on the fiber surfaces at different linear strains. Snapshots at linear strain of (a) 0.22, (b) 0.82, (c) 1.71, (d) 3.06, and (e) 6.39, respectively.

Corresponding to Fig. 10, in Fig. 11, the trace of the viscoelastic stress tensor is shown normalized by $\bar{τ}$ ⁠, where $\bar{τ}$ 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.

FIG. 11.

View large Download slide

Normalized trace of the viscoelastic stress tensor in the matrix for different linear strains. The snapshots are 2 D slices in the xz plane. For image interpretation and clarity, fibers close to the sliced plane are also shown. Snapshots at linear strain of (a) 0.22, (b) 0.82, (c) 1.71, (d) 3.06, and (e) 6.39, respectively.

C. Fiber orientation kinetics

In this section, the fiber orientation kinetics are investigated in uniaxial extension as a function of the Weissenberg number, strain hardening, and fiber volume fraction for suspensions with a viscoelastic matrix. The matrix, either LDPE or LLDPE, is modeled using a multi-mode Giesekus model, as described in Sec. VI A. RVE's containing 16 fibers are considered, since this optimizes the meaningful averages of the fiber orientation and the bulk stress as a function of computational time and number of realizations per parameter set, see the Appendix. The fibers in the RVE have an aspect ratio (L/D) of 10 with random initial positions and orientations. Hencky strain rates between 0.1 and 10 s⁻¹ are considered, and fiber volume fractions are varied between 1% and 10%. The initial domain size in these simulations is fixed, according to

L_{0, x} = L_{0, y} = L_{0, z} = \sqrt[3]{\frac{V_{p} N}{ϕ}},

(59)

where

V_{p}

is the volume of a single particle. The average anisotropy as a function of strain per parameter is shown by the solid lines in Fig. 12. The markers with errorbar in Fig. 12 represent the experimentally obtained anisotropy as a function of time, obtained in a previous work.¹ 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.¹

FIG. 12.

$Anisotropy of the fibers in the bulk as a function of linear strain. The solid lines are simulation results, and the markers with errorbar are experimental results obtained by Egelmeers et al.1 (a) Influence of the Weissenberg number (Hencky strain rate). The considered composites have an LDPE matrix and a fiber volume fraction of 1%, (b) influence of strain hardening. A Hencky strain rate of 1 s−1 is applied, the fiber volume fraction is 1%, and the matrix is either LDPE or LLDPE, and (c) influence of fiber volume fraction. Here, a Hencky strain rate of 1 s−1 is applied for an LDPE matrix.$

View large Download slide

Anisotropy of the fibers in the bulk as a function of linear strain. The solid lines are simulation results, and the markers with errorbar are experimental results obtained by Egelmeers et al.¹ (a) Influence of the Weissenberg number (Hencky strain rate). The considered composites have an LDPE matrix and a fiber volume fraction of 1%, (b) influence of strain hardening. A Hencky strain rate of 1 s⁻¹ is applied, the fiber volume fraction is 1%, and the matrix is either LDPE or LLDPE, and (c) influence of fiber volume fraction. Here, a Hencky strain rate of 1 s⁻¹ is applied for an LDPE matrix.

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⁻¹. 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⁻¹ 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⁻¹ 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. 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 ψ₀ and fixed θ₀ 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⁻¹ 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 θ₀ 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 $θ = cos (1 / \sqrt{3}) \approx 54.7$ °. This is clearly visible as a dip in viscosity for fibers that start with $θ_{0} > 54.7$ °. Up to 1 s, the strain hardening in FEM simulations is purely a result of fiber orientation kinetics.

FIG. 13.

$Transient uniaxial extensional viscosity for different initial fiber orientation angles θ 0 = arccos ( p z ). The solid lines correspond to Newtonian direct numerical FEM simulations. The dotted lines in (a) correspond to the viscosity prediction for the Dinh and Armstrong model with N p ϕ = 1.93. The dashed lines in (b) correspond to the viscosity computed using the constitutive model of Eq. (55) with N i = 2 , N p = 35, and N s = 0. (c) Normalized transient uniaxial extensional viscosity from FEM simulations with different numbers of fibers per simulation, and fixed initial fiber orientation angle θ 0 = 0°. All simulations are performed with a fiber volume fraction of 10% and fiber aspect ratio of 10.$

View large Download slide

Transient uniaxial extensional viscosity for different initial fiber orientation angles $θ_{0} = arccos (p_{z})$ ⁠. The solid lines correspond to Newtonian direct numerical FEM simulations. The dotted lines in (a) correspond to the viscosity prediction for the Dinh and Armstrong model with $N_{p} ϕ = 1.93$ ⁠. The dashed lines in (b) correspond to the viscosity computed using the constitutive model of Eq. (55) with $N_{i} = 2, N_{p} = 35$ ⁠, and $N_{s} = 0$ ⁠. (c) Normalized transient uniaxial extensional viscosity from FEM simulations with different numbers of fibers per simulation, and fixed initial fiber orientation angle $θ_{0} = 0$ °. All simulations are performed with a fiber volume fraction of 10% and fiber aspect ratio of 10.

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}, N_{p}$ ⁠, and $N_{s}$ ⁠) have been derived by Dinh and Armstrong²⁶ 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 (⁠ $ϕ$ ⁠) of 10%, $N_{p} ϕ = 1.93$ ⁠, and $N_{s} ϕ ≪ N_{p} ϕ$ ⁠, 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 ψ₀ and a fixed θ₀ 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 $θ_{0} = 60$ ° and $θ_{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 (⁠ $θ_{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, $ϕ$ ⁠, 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, N_{p} = {(L / D)}^{1.55}$ ⁠, and $N_{s} = 0$ (dashed lines), implying that up to the concentrated regime, $ϕ \approx 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.

FIG. 14.

$Transient uniaxial extensional viscosity for different initial fiber orientation angles θ 0 = arccos ( p z ). The solid lines correspond to Newtonian direct numerical FEM simulations. The dashed lines correspond to the constitutive model of Eq. (55), with N i = 2 , N p = ( L / D ) 1.55, and N s = 0. (a) Fiber volume fraction of 7% and fiber aspect ratio of 10, (b) fiber volume fraction of 4% and fiber aspect ratio of 10, (c) fiber volume fraction of 2% and fiber aspect ratio of 10, (d) fiber volume fraction of 10% and fiber aspect ratio of 7, (e) fiber volume fraction of 10% and fiber aspect ratio of 4, and (f) fiber volume fraction of 10% and fiber aspect ratio of 2.$

View large Download slide

Transient uniaxial extensional viscosity for different initial fiber orientation angles $θ_{0} = arccos (p_{z})$ ⁠. The solid lines correspond to Newtonian direct numerical FEM simulations. The dashed lines correspond to the constitutive model of Eq. (55), with $N_{i} = 2, N_{p} = {(L / D)}^{1.55}$ ⁠, and $N_{s} = 0$ ⁠. (a) Fiber volume fraction of 7% and fiber aspect ratio of 10, (b) fiber volume fraction of 4% and fiber aspect ratio of 10, (c) fiber volume fraction of 2% and fiber aspect ratio of 10, (d) fiber volume fraction of 10% and fiber aspect ratio of 7, (e) fiber volume fraction of 10% and fiber aspect ratio of 4, and (f) fiber volume fraction of 10% and fiber aspect ratio of 2.

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.¹ 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⁻¹, 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⁻¹ 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.

FIG. 15.

$Transient uniaxial extensional viscosity as a function of time for LDPE suspensions. The solid lines are simulation results, and the markers with error bars are experimental results obtained by Egelmeers et al.1 The dashed line indicates the region where the numerical results are affected by the time step choice, (a) rheological properties for LDPE with a fiber volume fraction of 1% at different Hencky strain rates, and (b) rheological properties for LDPE at a Hencky strain rate of 1 s−1 at different fiber volume fractions.$

View large Download slide

Transient uniaxial extensional viscosity as a function of time for LDPE suspensions. The solid lines are simulation results, and the markers with error bars are experimental results obtained by Egelmeers et al.¹ The dashed line indicates the region where the numerical results are affected by the time step choice, (a) rheological properties for LDPE with a fiber volume fraction of 1% at different Hencky strain rates, and (b) rheological properties for LDPE at a Hencky strain rate of 1 s⁻¹ at different fiber volume fractions.

FIG. 16.

View large Download slide

Transient uniaxial extensional viscosity as a function of time for LLDPE suspensions. The solid lines are simulation results, and the markers with error bar are experimental results obtained by Egelmeers et al.¹ The considered Hencky strain rate in all these simulations is 1 s⁻¹.

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⁻¹ 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 - 53$ ⁠), strain hardening, and fiber volume fraction (⁠ $ϕ = 1 % - 10 %$ ⁠), which agrees very well with experimental results obtained in our previous study.¹ 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,⁵⁵ who investigated orientation kinetics in planar extension far in the concentrated regime (⁠ $ϕ ≫ D / L$ ⁠). High agreement was also obtained between the rheology predicted by the numerical framework and the experimental results obtained in our previous study,¹ 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.⁵ 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 (⁠ $ϕ \approx 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 $η_{p} = 10^{5}$ Pa s, λ = 1 s, $α = 0.005, η_{sol} = 0$ Pa s is used, and a Hencky strain rate of ${\dot{ε}}_{H} = 1$ s⁻¹ 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.

FIG. 17.

View large Download slide

RVE statistics for different amounts of particles in the domain. The solid lines are the average and the dashed lines are the standard deviations. (a) Anisotropy of the fibers in the bulk as a function of linear strain, and (b) transient uniaxial extensional viscosity as a function of time.

REFERENCES

1.

T. R. N.

Egelmeers

,

N. O.

Jaensson

,

P. D.

Anderson

, and

R. M.

Cardinaels

, “

In-situ experimental investigation of fiber orientation kinetics during uniaxial extensional flow of polymer composites

,”

J. Rheol.

68

,

171–185

(

2024

).

https://doi.org/10.1122/8.0000749

Google Scholar

Crossref

2.

E. J.

Hinch

and

L. G.

Leal

, “

Time-dependent shear flows of a suspension of particles with weak Brownian rotations

,”

J. Fluid Mech.

57

(

4

),

753

–

767

(

1973

).

https://doi.org/10.1017/S0022112073001990

Google Scholar

Crossref

3.

S.

Yashiro

,

T.

Okabe

, and

K.

Matsushima

, “

A numerical approach for injection molding of short-fiber-reinforced plastics using a particle method

,”

Adv. Compos. Mater.

20

(

6

),

503

–

517

(

2011

).

https://doi.org/10.1163/092430411X584423

Google Scholar

Crossref

4.

M.

Yamanoi

,

J.

Maia

, and

T.

Kwak

, “

Analysis of rheological properties of fibre suspensions in a Newtonian fluid by direct fibre simulation. Part 2: Flexible fibre suspensions

,”

J. Non-Newtonian Fluid Mech.

165

(

19–20

),

1064

–

1071

(

2010

).

https://doi.org/10.1016/j.jnnfm.2010.05.004

Google Scholar

Crossref

5.

M.

Yamanoi

and

J.

Maia

, “

Analysis of rheological properties of fiber suspensions in a Newtonian fluid by direct fiber simulation. Part 3: Behavior in uniaxial extensional flows

,”

J. Non-Newtonian Fluid Mech.

165

(

23–24

),

1682

–

1687

(

2010

).

https://doi.org/10.1016/j.jnnfm.2010.09.006

Google Scholar

Crossref

6.

M.

Perez

,

A.

Scheuer

,

E.

Abisset-Chavanne

,

F.

Chinesta

, and

R.

Keunings

, “

A multi-scale description of orientation in simple shear flows of confined rod suspensions

,”

J. Non-Newtonian Fluid Mech.

233

,

61

–

74

(

2016

).

https://doi.org/10.1016/j.jnnfm.2016.01.011

Google Scholar

Crossref

7.

M.

Do-Quang

,

G.

Amberg

,

G.

Brethouwer

, and

A.

Johansson

, “

Simulation of finite-size fibers in turbulent channel flows

,”

Phys. Rev. E

89

(

1

),

013006

(

2014

).

https://doi.org/10.1103/PhysRevE.89.013006

Google Scholar

Crossref

8.

G.

D'Avino

,

M. A.

Hulsen

,

F.

Greco

, and

P. L.

Maffettone

, “

Bistability and metabistability scenario in the dynamics of an ellipsoidal particle in a sheared viscoelastic fluid

,”

Phys. Rev. E

89

(

4

),

043006

(

2014

).

https://doi.org/10.1103/PhysRevE.89.043006

Google Scholar

Crossref

9.

N. O.

Jaensson

,

M. A.

Hulsen

, and

P. D.

Anderson

, “

Simulations of the start-up of shear flow of 2d particle suspensions in viscoelastic fluids: Structure formation and rheology

,”

J. Non-Newtonian Fluid Mech.

225

,

70

–

85

(

2015

).

https://doi.org/10.1016/j.jnnfm.2015.09.006

Google Scholar

Crossref

10.

S. G.

Advani

and

C. L.

Tucker

, “

The use of tensors to describe and predict fiber orientation in short fiber composites

,”

J. Rheol.

31

,

751–784

(

1987

).

https://doi.org/10.1122/1.549945

Google Scholar

Crossref

11.

J. H.

Phelps

and

C. L.

Tucker

, “

An anisotropic rotary diffusion model for fiber orientation in short-and long-fiber thermoplastics

,”

J. Non-Newtonian Fluid Mech.

156

(

3

),

165

–

176

(

2009

).

https://doi.org/10.1016/j.jnnfm.2008.08.002

Google Scholar

Crossref

12.

J.

Wang

,

J. F.

O'Gara

, and

C. L.

Tucker

, “

An objective model for slow orientation kinetics in concentrated fiber suspensions: Theory and rheological evidence

,”

J. Rheol.

52

(

5

),

1179

–

1200

(

2008

).

https://doi.org/10.1122/1.2946437

Google Scholar

Crossref

13.

H.

Chen

,

P.

Wapperom

, and

D. G.

Baird

, “

The use of flow type dependent strain reduction factor to improve fiber orientation predictions for an injection molded center-gated disk

,”

Phys. Fluids

31

(

12

),

123105

(

2019

).

https://doi.org/10.1063/1.5129679

Google Scholar

Crossref

14.

Z.

Ouyang

,

E.

Bertevas

,

L.

Parc

,

B. C.

Khoo

,

N.

Phan-Thien

,

J.

Férec

, and

G.

Ausias

, “

A smoothed particle hydrodynamics simulation of fiber-filled composites in a non-isothermal three-dimensional printing process

,”

Phys. Fluids

31

(

12

),

123102

(

2019

).

https://doi.org/10.1063/1.5130711

Google Scholar

Crossref

15.

J.

Férec

,

E.

Bertevas

,

B. C.

Khoo

,

G.

Ausias

, and

N.

Phan-Thien

, “

Rigid fiber motion in slightly non-Newtonian viscoelastic fluids

,”

Phys. Fluids

33

(

10

),

103320

(

2021

).

https://doi.org/10.1063/5.0064191

Google Scholar

Crossref

16.

G. B.

Jeffery

, “

The motion of ellipsoidal particles immersed in a viscous fluid

,”

Proc. R. Soc. A

102

(

715

),

161

–

179

(

1922

).

https://doi.org/10.1098/rspa.1922.0078

Google Scholar

Crossref

17.

C. L.

Tucker

,

Fundamentals of Fiber Orientation: Description, Measurement and Prediction

(

Carl Hanser Verlag GmbH & Co. KG

,

2022

).

Google Scholar

Crossref

18.

S. K.

Kugler

,

A.

Kech

,

C.

Cruz

, and

T.

Osswald

, “

Fiber orientation predictions-A review of existing models

,”

J. Compos. Sci.

4

(

2

),

69

(

2020

).

https://doi.org/10.3390/jcs4020069

Google Scholar

Crossref

19.

A.

El-Rahman

and

C. L.

Tucker

, “

Mechanics of random discontinuous long-fiber thermoplastics. Part II: Direct simulation of uniaxial compression

,”

J. Rheol.

57

(

5

),

1463

–

1489

(

2013

).

https://doi.org/10.1122/1.4818804

Google Scholar

Crossref

20.

G. G.

Lipscomb

,

M. M.

Denn

,

D. U.

Hur

, and

D. V.

Boger

, “

The flow of fiber suspensions in complex geometries

,”

J. Non-Newtonian Fluid Mech.

26

(

3

),

297

–

325

(

1988

).

https://doi.org/10.1016/0377-0257(88)80023-5

Google Scholar

Crossref

21.

B. E.

VerWeyst

and

C. L.

Tucker

, “

Fiber suspensions in complex geometries: Flow/orientation coupling

,”

Can. J. Chem. Eng.

80

(

6

),

1093

–

1106

(

2002

).

https://doi.org/10.1002/cjce.5450800611

Google Scholar

Crossref

22.

R.

Tanner

, “

Aspects of non-colloidal suspension rheology

,”

Phys. Fluids

30

(

10

),

101301

(

2018

).

https://doi.org/10.1063/1.5047535

Google Scholar

Crossref

23.

R.

Tanner

, “

Rheology of noncolloidal suspensions with non-Newtonian matrices

,”

J. Rheol.

63

(

4

),

705

–

717

(

2019

).

https://doi.org/10.1122/1.5085363

Google Scholar

Crossref

24.

G. K.

Batchelor

, “

The stress generated in a non-dilute suspension of elongated particles by pure straining motion

,”

J. Fluid Mech.

46

(

4

),

813

–

829

(

1971

).

https://doi.org/10.1017/S0022112071000879

Google Scholar

Crossref

25.

C. J. S.

Petrie

, “

The rheology of fibre suspensions

,”

J. Non-Newtonian Fluid Mech.

87

(

2–3

),

369

–

402

(

1999

).

https://doi.org/10.1016/S0377-0257(99)00069-5

Google Scholar

Crossref

26.

S. M.

Dinh

and

R. C.

Armstrong

, “

A rheological equation of state for semiconcentrated fiber suspensions

,”

J. Rheol.

28

(

3

),

207

–

227

(

1984

).

https://doi.org/10.1122/1.549748

Google Scholar

Crossref

27.

R. B.

Pipes

,

J. W. S.

Hearle

,

A. J.

Beaussart

,

A. M.

Sastry

, and

R. K.

Okine

, “

A constitutive relation for the viscous flow of an oriented fiber assembly

,”

J. Compos. Mater.

25

(

9

),

1204

–

1217

(

1991

).

https://doi.org/10.1177/002199839102500907

Google Scholar

Crossref

28.

R. B.

Pipes

,

D. W.

Coffin

,

S. F.

Shuler

, and

P.

Šimáček

, “

Non-Newtonian constitutive relationships for hyperconcentrated fiber suspensions

,”

J. Compos. Mater.

28

(

4

),

343

–

351

(

1994

).

https://doi.org/10.1177/002199839402800404

Google Scholar

Crossref

29.

B.

Souloumiac

and

M.

Vincent

, “

Steady shear viscosity of short fibre suspensions in thermoplastics

,”

Rheol. Acta

37

,

289

–

298

(

1998

).

https://doi.org/10.1007/s003970050116

Google Scholar

Crossref

30.

J.

Férec

,

M. C.

Heuzey

,

J.

Pérez-González

,

L.

Vargas

,

G.

Ausias

, and

P. J.

Carreau

, “

Investigation of the rheological properties of short glass fiber-filled polypropylene in extensional flow

,”

Rheol. Acta

48

,

59

–

72

(

2009

).

https://doi.org/10.1007/s00397-008-0309-9

Google Scholar

Crossref

31.

G.

Ausias

,

X. J.

Fan

, and

R. I.

Tanner

, “

Direct simulation for concentrated fibre suspensions in transient and steady state shear flows

,”

J. Non-Newtonian Fluid Mech.

135

(

1

),

46

–

57

(

2006

).

https://doi.org/10.1016/j.jnnfm.2005.12.009

Google Scholar

Crossref

32.

A. H.

Wagner

,

D. M.

Kalyon

,

R.

Yazici

, and

T. J.

Fiske

, “

Uniaxial extensional flow behavior of a glass fiber-filled engineering plastic

,”

J. Reinf. Plast. Compos.

22

(

4

),

327

–

337

(

2003

).

https://doi.org/10.1177/0731684403022004280

Google Scholar

Crossref

33.

J.

Wang

,

W.

Yu

,

C.

Zhou

,

Y.

Guo

,

W.

Zoetelief

, and

P.

Steeman

, “

Elongational rheology of glass fiber-filled polymer composites

,”

Rheol. Acta

55

,

833

–

845

(

2016

).

https://doi.org/10.1007/s00397-016-0960-5

Google Scholar

Crossref

34.

C. W.

Macosko

,

Rheology: Principles, Measurements, and Applications

(

Wiley-VCH

,

1994

).

Google Scholar

35.

H.

Giesekus

, “

A simple constitutive equation for polymer fluids based on the concept of deformation-dependent tensorial mobility

,”

J. Non-Newtonian Fluid Mech.

11

(

1–2

),

69

–

109

(

1982

).

https://doi.org/10.1016/0377-0257(82)85016-7

Google Scholar

Crossref

36.

W. R.

Hwang

,

M. A.

Hulsen

, and

H. E. H.

Meijer

, “

Direct simulation of particle suspensions in sliding bi-periodic frames

,”

J. Comput. Phys.

194

(

2

),

742

–

772

(

2004

).

https://doi.org/10.1016/j.jcp.2003.09.023

Google Scholar

Crossref

37.

H. H.

Hu

,

N. A.

Patankar

, and

M.

Zhu

, “

Direct numerical simulations of fluid–solid systems using the arbitrary Lagrangian–Eulerian technique

,”

J. Comput. Phys.

169

(

2

),

427

–

462

(

2001

).

https://doi.org/10.1006/jcph.2000.6592

Google Scholar

Crossref

38.

G. K.

Batchelor

, “

The stress system in a suspension of force-free particles

,”

J. Fluid Mech.

41

(

3

),

545

–

570

(

1970

).

https://doi.org/10.1017/S0022112070000745

Google Scholar

Crossref

39.

W. R.

Hwang

,

M. A.

Hulsen

, and

H. E. H.

Meijer

, “

Direct simulations of particle suspensions in a viscoelastic fluid in sliding bi-periodic frames

,”

J. Non-Newtonian Fluid Mech.

121

(

1

),

15

–

33

(

2004

).

https://doi.org/10.1016/j.jnnfm.2004.03.008

Google Scholar

Crossref

40.

R.

Guénette

and

M.

Fortin

, “

A new mixed finite element method for computing viscoelastic flows

,”

J. Non-Newtonian Fluid Mech.

60

(

1

),

27

–

52

(

1995

).

https://doi.org/10.1016/0377-0257(95)01372-3

Google Scholar

Crossref

41.

A. C. B.

Bogaerds

,

A. M.

Grillet

,

G. W. M.

Peters

, and

F. P. T.

Baaijens

, “

Stability analysis of polymer shear flows using the extended pom-pom constitutive equations

,”

J. Non-Newtonian Fluid Mech.

108

(

1

),

187

–

208

(

2002

).

https://doi.org/10.1016/S0377-0257(02)00130-1

Google Scholar

Crossref

42.

A. N.

Brooks

and

T. J. R.

Hughes

, “

Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations

,”

Comput. Methods Appl. Mech. Eng.

32

(

1

),

199

–

259

(

1982

).

https://doi.org/10.1016/0045-7825(82)90071-8

Google Scholar

Crossref

43.

M. A.

Hulsen

,

R.

Fattal

, and

R.

Kupferman

, “

Flow of viscoelastic fluids past a cylinder at high Weissenberg number: Stabilized simulations using matrix logarithms

,”

J. Non-Newtonian Fluid Mech.

127

(

1

),

27

–

39

(

2005

).

https://doi.org/10.1016/j.jnnfm.2005.01.002

Google Scholar

Crossref

44.

R.

Fattal

and

R.

Kupferman

, “

Constitutive laws for the matrix-logarithm of the conformation tensor

,”

J. Non-Newtonian Fluid Mech.

123

(

2

),

281

–

285

(

2004

).

https://doi.org/10.1016/j.jnnfm.2004.08.008

Google Scholar

Crossref

45.

G.

Geuzaine

and

J. F.

Remacle

, “

Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities

,” Int. J.

Numer. Methods Eng.

79

(

11

),

1309

–

1331

(

2009

).

https://doi.org/10.1002/nme.2579

Google Scholar

Crossref

46.

O.

Schenk

,

K.

Gärtner

, and

W.

Fichtner

, “

Efficient sparse LU factorization with left-right looking strategy on shared memory multiprocessors

,”

BIT Numer. Math.

40

,

158

–

176

(

2000

).

https://doi.org/10.1023/A:1022326604210

Google Scholar

Crossref

47.

O.

Schenk

and

K.

Gärtner

, “

Solving unsymmetric sparse systems of linear equations with PARDISO

,”

Future Gener. Comput. Syst.

20

(

3

),

475

–

487

(

2004

).

https://doi.org/10.1016/j.future.2003.07.011

Google Scholar

Crossref

48.

O.

Schenk

and

K.

Gärtner

, “

On fast factorization pivoting methods for sparse symmetric indefinite systems

,”

Electron. Trans. Numer. Anal.

23

(

1

),

158

–

179

(

2006

).

Google Scholar

49.

G.

D'Avino

and

M. A.

Hulsen

, “

Decoupled second-order transient schemes for the flow of viscoelastic fluids without a viscous solvent contribution

,”

J. Non-Newtonian Fluid Mech.

165

(

23–24

),

1602

–

1612

(

2010

).

https://doi.org/10.1016/j.jnnfm.2010.08.007

Google Scholar

Crossref

50.

D. C.

Rapaport

,

The Art of Molecular Dynamics Simulation

(

Cambridge University Press

,

2004

).

Google Scholar

Crossref

51.

J.

Diebel

, “

Representing attitude: Euler angles, unit quaternions, and rotation vectors

,”

Matrix

58

(

15–16

),

1

–

35

(

2006

).

Google Scholar

52.

N. O.

Jaensson

,

M. A.

Hulsen

, and

P. D.

Anderson

, “

Stokes–Cahn–Hilliard formulations and simulations of two-phase flows with suspended rigid particles

,”

Comput. Fluids

111

,

1

–

17

(

2015

).

https://doi.org/10.1016/j.compfluid.2014.12.023

Google Scholar

Crossref

53.

R.

Takserman-Krozer

and

A.

Ziabicki

, “

Behavior of polymer solutions in a velocity field with parallel gradient. I. Orientation of rigid ellipsoids in a dilute solution

,”

J. Polym. Sci., Part A

1

(

1

),

491

–

506

(

1963

).

https://doi.org/10.1002/pol.1963.100010143

Google Scholar

Crossref

54.

C. L.

Tucker

, “

Flow regimes for fiber suspensions in narrow gaps

,”

J. Non-Newtonian Fluid Mech.

39

(

3

),

239

–

268

(

1991

).

https://doi.org/10.1016/0377-0257(91)80017-E

Google Scholar

Crossref

55.

G. M.

Lambert

and

D. G.

Baird

, “

Evaluating rigid and semiflexible fiber orientation evolution models in simple flows

,”

J. Manuf. Sci. Eng.

139

(

3

),

031012

(

2017

).

https://doi.org/10.1115/1.4034664

Google Scholar

Crossref

Numerical simulation of fiber orientation kinetics and rheology of fiber-filled polymers in uniaxial extension Open Access

I. INTRODUCTION

II. PROBLEM DESCRIPTION

A. Particles in the triperiodic domain

B. Balance equations

C. Boundary and initial conditions

D. Arbitrary Lagrange–Euler (ALE) formulation

E. Bulk stresses

F. Anisotropy

G. RVE statistics

III. NUMERICAL METHOD

A. Weak formulation

B. Spatial discretization

C. Time discretization

D. Remeshing and projection

IV. NUMERICAL MODEL VERIFICATION

A. Convergence study

B. Determination of the mesh size

V. CONSTITUTIVE MODEL FOR RHEOLOGY OF FIBER SUSPENSIONS

VI. RESULTS

A. Matrix constitutive model parameters

B. In situ deformation and stress fields

C. Fiber orientation kinetics

D. Bulk rheology

1. Newtonian matrix

2. Viscoelastic matrix

VII. CONCLUSION AND DISCUSSION

SUPPLEMENTARY MATERIAL

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

Author Contributions

DATA AVAILABILITY

APPENDIX: DETERMINATION OF THE RELATIVE RVE SIZE

REFERENCES

Supplementary Material

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

Numerical simulation of fiber orientation kinetics and rheology of fiber-filled polymers in uniaxial extension