In this study, we experimentally investigate the fiber orientation kinetics and rheology of fiber-filled polymer melts in shear flow. A novel setup is designed with custom-built bottom and top geometries that are mounted on a conventional rotational rheometer. Shear flow between parallel sliding plates is applied by vertical movement of the top geometry. The axial force measurement data of the rotational rheometer are used to determine the shear stress growth coefficient. The fiber orientation kinetics are measured in situ with this setup using small angle light scattering. We consider a non-Brownian experimental system with short glass fibers for the suspended phase ( L / D = 8 15) and different polyethylene based materials for the matrix phase. The fiber orientation kinetics are investigated as a function of fiber volume fraction ( ϕ = 1 %, 5%, and 10%) and as a function of the shear rate ( γ ˙ = 0.03, 0.55, and 5 s 1). Within the studied range, these parameters do not influence the fiber orientation kinetics, and a multiparticle model, based on Jeffery’s equation for single particles, can describe these kinetics. Our results show that, up to the concentrated regime ( ϕ D / L), fiber-fiber interactions do not influence the fiber orientation in shear flow. Finally, we investigate the shear stress growth coefficient of these composites and demonstrate that a simple rheological model for fiber composites, which assumes a constant, isotropic orientation distribution of the fibers, is able to describe the shear stress growth coefficient of the short fiber composite samples.

Fiber-filled polymer materials are increasingly used in the high-tech industry because of their unique combined properties, which, among others, are a high strength, a low weight, and a low manufacturing cost. The properties of processed fiber-filled products can be highly anisotropic, depending on the local fiber orientation. Therefore, it is of interest to develop an understanding of the relation between polymer flow during processing and fiber orientation. In our previous work [1,2], the fiber orientation kinetics in uniaxial extensional flow were analyzed. In the current work, the fiber orientation kinetics in shear flow will be investigated.

Particle orientation kinetics in shear flow were first investigated by Jeffery in 1922 [3]. He developed a theoretical framework describing the orientation evolution of a single ellipsoidal particle, wherein he assumed a Newtonian fluid and a Stokes flow field. This theory can be applied in the dilute regime, where no fiber-fiber interactions are present. Folgar and Tucker [4] were the first to develop a theoretical framework to describe the orientation kinetics for a group of fibers. The basis of the model is the theory as described by Jeffery but extended with phenomenological terms. The standard Folgar–Tucker model includes a phenomenological Brownian diffusion term to include fiber-fiber interactions that hinder periodic particle rotation [5,6]. This standard model is commonly extended with another phenomenological term, which is called the “slip factor,” and is required to better reflect the slower orientation kinetics experimentally observed for shear flows, when comparing to the standard Folgar–Tucker model [6–10]. The physical mechanism behind this slow orientation kinetics is currently unknown [11].

The most common experimental approaches to investigate the orientation kinetics of anisotropic particles (such as fibers) in shear flow include x-ray tomography [12], scattering techniques [7,13–19], microscopy [6,10,20–23], and rheology [8,9,24,25]. In the rheology approach, shear stress growth coefficients (in the remainder of this work, the term transient shear viscosity will be used instead of shear stress growth coefficient) and first normal stress growth functions of fiber suspensions are measured. In these works, first a negative preshear is applied to align fibers to some extent in the velocity direction, when subsequently shearing in the opposite direction the fibers tumble and this results in additional contributions to the viscosity and normal force, which can be used to determine fiber orientation kinetics. This is an indirect in situ method, which does not provide information about the orientation of the individual fibers, and is limited to concentrated fiber suspensions. Scattering techniques are, in general, also indirect in situ methods and allow to determine the fiber orientation kinetics starting from an isotropic sample with high acquisition rates, but this method only provides a 2D-averaged projection of the structure. The resulting anisotropy can, therefore, be interpreted in multiple manners. Microscopy methods do provide exact information about the individual fibers but are limited to dilute suspensions [6,20]. X-ray tomography can provide 3D information on individual fibers in a concentrated suspension, but capturing fast kinetics requires synchrotron radiation [12].

Extensive research has been conducted to investigate the orientation of anisotropic particles for suspensions in shear flow with low viscous, often Newtonian, matrices [17–19]. Interesting phenomena like log rolling and tumbling have been found at large strains. For fibers suspended in high viscous or viscoelastic matrices, various experimental characterization methods led to the following qualitatively similar conclusions. Orientation kinetics slow down with increasing fiber volume fraction and fiber aspect ratio [6–10]. Furthermore, for non-Brownian suspensions and for polypropylene and polybutene matrices, the kinetics are controlled by the strain, independent of the shear rate, within the shear rate ranges considered in these works [8,10]. Moreover, already for semiconcentrated fiber suspensions, the orientation kinetics are much slower than predicted by Jeffery’s equation and are, for concentrated suspensions, still evolving after 100 shear strain units [21,23,25]. Quantitatively, however, the orientation kinetics are not well understood, because for each experimental system (fiber aspect ratio, fiber flexibility, fiber volume fraction, matrix) different phenomenological parameters in the Folgar–Tucker model are required to describe the orientation kinetics.

The orientation of fibers also influences the rheology of fiber composites, which are generally measured using either rotational rheometers [21,26–30] or sliding plate rheometers [10,22,23,31–33]. In the case of rotational shear devices (cone-plate and parallel plate geometries), the presence of the curvilinear streamlines and the rigid fibers that prefer to align in the flow direction can lead to inaccurate rheological measurements [5,21]. Furthermore, the inhomogeneous shear field in the parallel plate geometries, in combination with the strain dependent fiber orientation kinetics [8,10], implies that the fiber orientation kinetics in the sample are a function of the radial coordinate, which can lead to additional fiber-fiber interactions and inaccurate rheological measurements. A cone and plate geometry does have a homogeneous shear field, but generally does not allow for gap control and the fibers in the samples will be preoriented after the composite being “squeezed” in the narrow gap during sample loading. Although Eberle et al. [21] used donut shaped samples to prevent squeezing effects at low radius, still, for the cone and plate geometry also boundary interactions or confinement effects can prevent fibers from rotating freely [34,35], which can influence the measurement. Confinement effects can also introduce shear banding phenomena, which further complicate the fiber orientation and rheology, as observed by Li et al. [36]. To investigate the transient shear rheology, sliding plate rheometers are preferred because they have a homogeneous shear field, and the gap can be controlled to prevent fiber-boundary interactions. The drawback of these devices is a limited maximum achievable strain.

Using a sliding plate rheometer, both Cieslinski et al. [23] and Ortman et al. [10,31] investigated the transient shear viscosity of polypropylene glass fiber suspensions. They both noticed that the zero shear viscosity is dependent on the initial fiber orientation. Furthermore, they found a fiber induced shear stress overshoot in the transient viscosity curve for a random initial orientation. To describe the transient viscosity using a rheological constitutive model, Ortman et al. [10] used an additional term in the expression for the stress tensor taking fiber flexibility into account, whereas Cieslinski et al. [23] used an additional term for fiber-fiber interaction. Both these works contributed greatly to understanding the rheological properties of long glass fiber-filled polymers. However, they did not consider the rheology of short glass fiber composites, for which the effect of orientation on the rheology may be smaller.

In the present study, the fiber orientation kinetics as a function of Weissenberg number, matrix rheology, and fiber volume fraction will be investigated for non-Brownian short glass fiber polymer composites, in situ, in the velocity–velocity gradient plane of a shear field. A novel setup is designed, which is based on a commercial rotational rheometer and applies shear between parallel sliding plates, and which measures the fiber orientation kinetics using small angle light scattering (SALS). The experimentally obtained fiber orientation kinetics will be compared to the fiber orientation kinetics obtained from a model based on Jeffery’s equation. Furthermore, the force resisting the movement of the parallel plate will be used to determine the viscosity of the fiber composites, which will be compared to the predictions of a simple constitutive equation.

In this section, first, a model is introduced that describes the anisotropy evolution of a polymer composite containing noninteracting rodlike particles in shear flow. This model is based on Jeffery’s equation [3], which is derived for Newtonian fluids. Afterward, a second model is presented, describing the transient rheology of fiber-filled composites.

The orientation of a single rodlike particle can be described using two angles: ψ p ( t ) and θ p ( t ), defined in Fig. 1. Figure 1 also shows that shear flow is applied in the x y-plane. The angle between the y-axis and the orientation of the particle projected in the velocity–velocity gradient plane of the shear field is the tumble angle, ψ p ( t ). The second angle, θ p ( t ), is the angle the particle makes with the z-axis. The orientation vector, p, can be written using these angles (after conversion of Cartesian coordinates to spherical coordinates),
(1)
with
(2)
(3)
(4)
Moreover, the velocity, velocity gradient, and vorticity directions are indicated with e x, e y, and e z. The rate of change of p of a particle in a Newtonian fluid is described by Jeffery’s equation and is given by [11]
(5)
where β is the particle shape factor, which is related to the geometry and aspect ratio of the rodlike particle. For a prolate ellipsoid, β = ( r e 2 1 ) / ( r e 2 + 1 ), with r e = L / D and L being the particle length and D being the particle diameter. For large aspect ratios, β approaches unity. Furthermore, D is the rate-of-deformation tensor and W is the vorticity tensor,
(6)
(7)
where u is the velocity vector. For simple shear, these tensors are W = 1 2 γ ˙ ( e x e y e y e x ) and D = 1 2 γ ˙ ( e y e x + e y e x ), where γ ˙ is the shear rate.
FIG. 1.

Schematic representation of the relevant angles to describe the orientation of a particle. The shear field is shown in the x y-plane, according to u = γ ˙ y e x.

FIG. 1.

Schematic representation of the relevant angles to describe the orientation of a particle. The shear field is shown in the x y-plane, according to u = γ ˙ y e x.

Close modal
After substitution of Eqs. (1)–(4) in Eq. (5), the resulting three equations can be solved for the unknown angles. After simplification, the differential equations for both angles are [11] as follows:
(8)
(9)
Note that ψ p = 90 ° in this study corresponds to ψ p = 0 ° in the book of Tucker [11], corresponding to perfect alignment of the fiber in the flow direction. The above equations can be integrated using the initial conditions ψ p ( t = 0 ) = ψ p , 0 and θ p ( t = 0 ) = θ p , 0 to yield
(10)
with C 1 being a constant,
(11)
(12)
where δ s = π if p y , 0 < 0 and δ s = 0 if p y , 0 0, and C 2 is an integration constant,
(13)

In the remainder of this work, only particles with an initial orientation in the velocity–velocity gradient plane ( p z , 0 = 0) are considered. For this case, θ p , 0 = 90 °, so that C 2 and as a result θ p ( t ) = 90 °. Thus, the particle orientation kinetics are fully determined by solving Eq. (10). Figure 2 shows the kinetics of the tumble angle ψ p for 20 particles with randomized initial orientation and r e = 10. The figure shows that tumbling occurs around ψ p = 0 °, 180 °, and 360 °, where particles align perpendicular to the velocity direction. The plateaus for ψ p = 90 °, 270 °, and 450 ° correspond to particles orienting perfectly in the velocity direction. Note that a tumbling angle ψ p ± k 180 °, with k an integer, describes the same particle orientation.

FIG. 2.

Tumble angle, ψ p, as a function of shear strain for different initial particle orientations, according to Jeffery’s equation [Eq. (10)]. The particles have an aspect ratio r e = 10. The inset indicates the shear plane and a particle oriented in the flow direction and a particle perpendicular to the flow direction.

FIG. 2.

Tumble angle, ψ p, as a function of shear strain for different initial particle orientations, according to Jeffery’s equation [Eq. (10)]. The particles have an aspect ratio r e = 10. The inset indicates the shear plane and a particle oriented in the flow direction and a particle perpendicular to the flow direction.

Close modal
In the multiparticle model, multiple particles are initially randomly oriented in the velocity–velocity gradient plane, and Eq. (10) is solved for each particle. This model, thus, neglects particle-particle interactions. Subsequently, at every time step the particle orientation tensor, P, is computed by averaging the dyadic product between the particle unit vectors, according to [4]
(14)
where N is the number of particles in the multiparticle model and p i is the orientation vector for particle i, without the p z component so that P is a 2 × 2 tensor. The eigenvalues and the corresponding eigenvectors of the particle orientation tensor are computed and are used to determine the anisotropy parameter, ξ m, and the average orientation angle, ψ m, as follows:
(15)
(16)
where ν m , max and ν m , min are the largest and smallest eigenvalues of the particle orientation tensor, respectively, and ι m , max is the eigenvector corresponding to the largest eigenvalue. Negative average orientation angles are shifted with 180 ° so that ψ m is in the range [ 0 , π].

In case all particles are aligned in the y-direction or the x-direction, the average orientation angle ψ m = 0 ° and ψ m = 90 °, respectively, while in both cases the anisotropy ξ m = 1. Initially, for a random particle orientation in the velocity–velocity gradient plane, ν m , max = ν m , min = 1 / 2, and the anisotropy parameter becomes ξ m = 0. The average orientation angle is, in the latter case, not well defined. Note, as opposed to elongational flow for which a similar model was derived in our previous work [1], in shear flow, the off-diagonal components of the particle orientation tensor are not always equal to zero and, therefore, the average orientation angle is relevant. By systematically increasing the number of particles in the multiparticle model, it was determined that 2000 particles are sufficient to obtain a reproducible average anisotropy and average orientation angle. For this amount of particles, the anisotropy and the average orientation angle as a function of strain are investigated for different particle aspect ratios, see Fig. 3. The figure shows that for all considered particle aspect ratios, the anisotropy first increases due to particles aligning in the velocity direction and subsequently decreases due to the tumbling motion of the individual particles. The higher the particle aspect ratio, the longer the particles remain aligned in the velocity direction, and the higher the anisotropy of the system. The figure also shows that for particle aspect ratios larger than 7, the initial orientation kinetics become independent of the particle aspect ratio. At these large aspect ratios, only the tumbling period is dependent on the particle aspect ratio. Although this model is derived for Newtonian fluids, the aim is to assess its applicability in describing particle orientation for viscoelastic fluids. Note, in case both initial particle orientation angles are random, the model prediction of the anisotropy in the velocity–velocity gradient plane increases slightly faster than when the particles are initially randomly oriented in the velocity–velocity gradient plane, whereas the average orientation angle in the velocity–velocity gradient plane is exactly identical for the two cases, see Sec. S1 of the supplementary material for a comparative figure.

FIG. 3.

(a) Anisotropy and (b) average orientation angle of particles with various aspect ratios, r e, as a function of shear strain, based on the multiparticle model [Eqs. (15) and (16)] with random initial orientation in the velocity–velocity gradient plane.

FIG. 3.

(a) Anisotropy and (b) average orientation angle of particles with various aspect ratios, r e, as a function of shear strain, based on the multiparticle model [Eqs. (15) and (16)] with random initial orientation in the velocity–velocity gradient plane.

Close modal
The starting point of the model describing the rheology of fiber composites is the general constitutive fiber model first proposed by Hinch and Leal [37] in 1973. Later, exact expressions were found to describe the parameters of this model by Dinh and Armstrong [38]. After neglecting the Brownian motion term in these models, a simplified expression for the transient shear viscosity is found, according to [11]
(17)
where ϕ is the fiber volume fraction, D x y = γ ˙ / 2 is the x y-component of the rate-of-deformation tensor and P x y x y is the component in the shear plane of the fourth order fiber orientation tensor, which is computed using a fiber orientation tensor P = 1 / 2 e x e x + 1 / 2 e y e y, complying with a random particle orientation distribution and the natural closure approximation. The experimental method, see Sec. III D, only allows for rheology measurements at low strains, where orientation effects are not significant yet. Therefore, in the model, P x y x y = 1 / 4 is constant. The isotropic viscosity, η I +, includes both the contribution of the matrix and that of the particles,
(18)
where η s , pure + is the transient shear viscosity of the unfilled matrix. The isotropic contribution of the particles to the viscosity is captured by parameter N i. In our previous work on fiber composites in uniaxial extension [2], we found by fitting a similar model on FEM simulations that N i = 2, and Hinch and Leal [39] also found N i = 2 for slender ellipsoids with an aspect ratio of r e 1. Lipscomb et al. [40] used the elliptic integrals of Giesekus [41] to exactly derive N i for ellipsoids with aspect ratios between 0 and infinity. They also found N i = 2 for large aspect ratios, and they showed that for spherical particles ( L / D = 1), the coefficient N i = 2.5, which is identical to the Einstein equation. The anisotropic contribution of the particles to the viscosity is captured with the particle number, N p, which can be written as [11,38,42]
(19)
Substitution of Eqs. (19) and (18) in Eq. (17) and simplification results into
(20)
In this model with a constant random particle orientation, the transient shear viscosity of fiber composites is, thus, only a function of the fiber volume fraction, the fiber aspect ratio, and the matrix rheology. In case particles are fully aligned in the flow direction or velocity gradient direction, the product P x y x y D x y = 0, independent of the closure approximation and, therefore, the anisotropic contribution of the fibers to the stress will be zero. In the remainder of this work, the fiber aspect ratio in the model, r e = 10, unless otherwise indicated.

In this section, first, the materials, and the sample preparation are explained. Then, the setup for measuring both the fiber orientation and the rheological properties during shear flow is discussed.

In this study, glass fibers (EPH80M-10A) are used as the suspended phase, which are kindly provided by Nippon Electric Glass. The fibers have a diameter of 11  μm and a volume-averaged length of 111  μm, which are determined on extruded and compression-molded fiber composite samples using x-ray tomography. The fiber length distribution in the prepared samples can be found in the supplementary material of our previous work [1]. Furthermore, the fibers are amino-silane-coated, see our previous study [1]. The matrix materials are low density polyethylene (LDPE), linear low density polyethylene (LLDPE), and maleic anhydride grafted low density polyethylene (LDPE-g-MAH, commercial name Yparex 0H31). LDPE ( M w = g mol 1, M n = g mol 1) has a viscosity-averaged relaxation time of 53 s at 150 °C and exhibits strain hardening in extension. LLDPE ( M w = g mol 1, M n = g mol 1) has a viscosity-averaged relaxation time of 0.62 s at 150  °C and does not exhibit strain hardening in extension. LDPE-g-MAH has a melt flow index of 6.0 g/10 min at 190  °C and 2.16 kg and is kindly provided by The Compound Company. The rheology of LDPE and LLDPE is fully characterized in our previous studies [1,2]. LDPE-g-MAH is expected to covalently bond with the amino-silane coating on the glass fibers for better matrix-fiber adhesion, which prevents slip. According to the manufacturer, the amount of maleic anhydride is around 2%–3% based on mass. For this experimental system, the ratio between convection effects and diffusion effects (or Péclet number) is of order 10 8, indicating that Brownian forces are negligible.

First, polyethylene and glass fibers are melt mixed with a DSM Micro 15cc twin screw compounder. The mixing is performed using a rotational speed of 50 rpm (8 min) at 230  °C. To reduce air bubbles in the extrudate, the rotational speed is lowered to 10 rpm after 7 min.

The extrudate is then cut into small granules and is compression molded at 160 °C to form thin sheets. To achieve maximal randomization of fiber orientation, these thin sheets are further cut into smaller pieces of approximately 2  × 2  mm 2 and compression molded once more into thin sheets. Subsequently, to prevent shear buckling in our experimental setup, the thickness of the samples is increased to 1.4 mm by compression molding pure matrix material on top of the thin sheet composites. This way, high local fiber volume fractions can be achieved while still enabling light scattering experiments. During this compression molding step, the smoothness of the bottom and top surfaces is increased by a protective Kapton film. The increased smoothness is required to perform scattering experiments at room temperature. Finally, the sheets are cut into samples with dimensions of approximately 20 × 8  mm 2 (length × width). During all compression molding cycles, the following pressure-time profile is used: First, no pressure is applied (3 min); second, a pressure of 20 kN is applied (3 min); and third, a pressure of 40 kN is applied (5 min).

Fiber volume fractions of 1%, 5%, and 10% are made (2.6%, 13%, and 26% fiber weight fraction, respectively). The fiber orientation is measured using small angle light scattering, which limits the amount of fibers in the sample. Therefore, the thickness of the sample part that is filled with fibers is a function of the fiber volume fraction. For fiber volume fractions of 1%, 5%, and 10%, the thickness of the sample filled with fibers is 1.4, 0.44, and 0.22 mm, respectively. For the rheological experiments of the composites, samples with a thickness of 1.4 mm are used. These samples are fully filled with fibers.

A Nanotom 160 NF (General Electric/Phoenix) with a voltage of 60 kV, a current of 310  μA, and a voxel size between 1.3 and 2.0  μ m 3 is used to perform CT scans. GeoDict [43] is used to analyze the scans. The software package GeoDict allows for the determination of the start and end coordinates of individual fibers in 3D. These coordinates are used to compute the fiber orientation vectors and subsequently the fiber orientation tensor, similar to Eq. (14). After verifying that the thickness components of the fiber orientation tensor are close to zero, see [1], this 3 × 3 tensor can be written as a 2 × 2 tensor, and the anisotropy can be computed using the eigenvalues of the 2 × 2 fiber orientation tensor,
(21)
where ν CT , min and ν CT , max are the smallest and largest eigenvalues of the 2 × 2 fiber orientation tensor resulting from the CT scan.

A custom setup is built in-house for the in situ determination of the fiber orientation kinetics in shear flow using SALS. In our study, the chosen setup is optimized to capture fiber orientation in the velocity–velocity gradient plane, which is essential for observing the orientation angle during tumbling. Commercial rheo-SALS setups either restrict SALS measurements to the velocity–vorticity plane (in a parallel plate setup) or are not suitable for highly viscous polymer melts (in a Couette setup). The setup is schematically shown in Fig. 4. From right to left, the setup consists of a Neon-Helium laser with a wavelength of 632.8 nm and a spot size of about 1  mm 2, an Anton Paar MCR 502 rheometer, a white detector plate, and a Hamamatsu C14400 camera having the capability of capturing 100 frames per second. The rheometer is equipped with a Convection Temperature control Device (CTD) 300/GL, which is a special oven with optically flat glass windows perpendicular to the laser source. A cuvette with quartz glass windows, shown in Fig. 5, is placed in the oven. The cuvette is equipped with an aluminum stationary pillar, serving as the fixed sample holder. The distance from the stationary pillar to the oven’s centerline is adjustable, allowing to vary the sample gap height. Moreover, the top geometry of the Anton Paar rheometer consists of a shaft used for disposable plate-plate rheology, which is extended with a custom shaped aluminum rectangular prism, serving as the moving sample holder. The cuvette is filled with silicon oil, having a viscosity of 10 mPa s. Immersing the samples in the oil is required to prevent unwanted scattering due to surface roughness generated when shearing the fiber-filled samples. The laser is hitting the sample in the x y-plane and in the y-direction, it is exactly in the center between both sliding plates, as shown in (9) in Fig. 5. The ratio between the gap size and the laser spot diameter is about 1.8; therefore, no boundary effects are expected to occur when capturing the orientation kinetics. The distance from the sample to the detector plate is approximately 32.5 cm so that characteristic distances between 3 and 9  μm are visible in the scattering pattern ( q-range from 7.0 × 10 4 to 2.1 × 10 3 nm 1). The considered q-range is just above q = 2 π / D. Whereas shape anisotropy can be detected for q < 2 π / D, our data cover a q-range where the structure and form factor decay and there is a transition toward the high q-range determined by the surface state of the fibers. Hence, although the fibers appear relatively smooth in scanning electron microscopy (SEM) images (Figs. S13 and S14 in the supplementary material), any fiber surface anisotropy will also contribute to the measured scattering signal. The considered q-range is, however, too limited to extract detailed structural information from the intensity versus q data.

FIG. 4.

Schematic depiction of the shear-SALS setup. The most relevant components of the setup, from right to left, are as follows: a Neon-Helium laser (1), a convection oven (2) with an oil-filled cuvette (3), a shaft for disposable geometries with a custom-made prism geometry (4), a detector plate (5), and a camera (6).

FIG. 4.

Schematic depiction of the shear-SALS setup. The most relevant components of the setup, from right to left, are as follows: a Neon-Helium laser (1), a convection oven (2) with an oil-filled cuvette (3), a shaft for disposable geometries with a custom-made prism geometry (4), a detector plate (5), and a camera (6).

Close modal
FIG. 5.

Schematic representation of the construction in the convection oven to apply shear on the sample, consisting of an oil-filled cuvette (3), a fixed pillar used to set the sample gap (7), the composite sample (8), a vertically moving pillar, which is a rheometer shaft with a rectangular prism attached to it (4). The laser location is indicated with (9). (a) The initial geometry of the sample and initial position of the moving pillar and (b) the geometry of the sample and the corresponding position of the moving pillar after applying a certain amount of shear strain. A more detailed visualization of the sheared samples is shown in Fig. 9.

FIG. 5.

Schematic representation of the construction in the convection oven to apply shear on the sample, consisting of an oil-filled cuvette (3), a fixed pillar used to set the sample gap (7), the composite sample (8), a vertically moving pillar, which is a rheometer shaft with a rectangular prism attached to it (4). The laser location is indicated with (9). (a) The initial geometry of the sample and initial position of the moving pillar and (b) the geometry of the sample and the corresponding position of the moving pillar after applying a certain amount of shear strain. A more detailed visualization of the sheared samples is shown in Fig. 9.

Close modal

In all these experiments, at least three repeat measurements are performed per parameter set. The maximum achievable strain in these experiments is about γ = 10. The fiber orientation kinetics for composites with either LDPE or LLDPE as matrix material are determined at a temperature of 150 °C, while for composites with LDPE-g-MAH as matrix material the temperature is set to 140 °C.

Using clamps, the composite sample is attached to both the stationary pillar and the rectangular prism of the top geometry, see Fig. 5. For clarity, the figure does not include these clamps. Shear flow is applied by downward vertical movement of the top geometry, which is also schematically shown in Fig. 5. The whole top geometry is free to rotate, leading to small variations in the distance between the stationary pillar and the rectangular prism, which determines the sample gap. Typical values of the sample gap range from 1.6 to 2.2 mm. Therefore, in each experiment, the sample gap is measured using a camera and subsequently, the downward velocity of the top geometry is adjusted so that the desired shear rate is applied on the sample. Furthermore, during the sample loading and the downward motion, the rotation speed of the top geometry is set to zero. The vertical acceleration of the top geometry is restricted to 10 mm  s 2, resulting in a gradual increase in strain rate as a function of time. This limitation is only significant for the highest strain rate. To mitigate this effect, a negative prestrain ( γ r = 2) is applied (only for the highest shear rate), displacing the top geometry vertically upward, aligning fibers to some extent. Finally, after relaxation of the stress, the regular downward velocity is applied onto the sample, ensuring that at the initial vertical position of the top geometry (with random fiber orientation), the desired high shear rate is applied.

The Weissenberg number in the matrix of the composites can be computed, according to
(22)
where λ c is the viscosity-averaged relaxation time of the matrix material, which was found to be 53 and 0.62 s for LDPE and LLDPE, respectively, at 150  °C, see our previous work [1]. The considered shear rates in the experiments in the current paper are γ ˙ = 0.03, 0.5, and 5.0  s 1, which then corresponds to a viscosity-averaged Weissenberg number ranging from Wi = 1.6 to 265 for LDPE and Wi = 0.02 to 3.1 for LLDPE. Based on the frequency sweep data of our previous work [1], the reptation and Rouse time of LDPE and LLDPE can be estimated, and accordingly the reptation and Rouse time based Weissenberg numbers, Wi rep and Wi rouse, see Sec. S2 of the supplementary material. For LDPE and LLDPE, the reptation time based Weissenberg number ranges from 8.1 to 1350 and from 0.51 to 85, respectively. For Wi rep > 1, polymer chain orientation is expected, which is the case for all shear rates considered in the current paper except the lowest shear rate for composites with an LLDPE matrix. For LDPE and LLDPE, the Rouse time based Weissenberg number ranges from 8.1 to 1350 ( O ( Wi rep )) and from 0.0011 to 0.18, respectively. For Wi rouse > 1, polymer chain stretching is expected, which is the case for all composites with an LDPE matrix.

1. Shear rheology in shear-SALS setup

The downward motion of the moving pillar imposes a shear flow on the sample. The sample’s resistance against this flow is determined by the viscosity, which is measured in the shear-SALS setup using the axial force sensor of the Anton Paar MCR 502 rheometer, according to
(23)
where d is the minimum distance between the stationary pillar and the rectangular prism attached to the shaft of the rheometer; u y is the downward velocity of the top geometry; A is the cross-sectional area of the sample (sample thickness × sample height), F is the axial force, measured during the downward motion; F ns is the axial force measured during the downward motion without attachment of a sample; and C γ ˙ is a correction factor on the applied shear rate, which will be further discussed in Secs. IV and V A 2. For this paper, the rheological measurements are performed in the cuvette without silicon oil, hence F ns 0, and in situ SALS experiments are performed separately with oil. The transient shear viscosity data are only accurate up to a strain of about γ = 2.6, as demonstrated in Sec. V B 2.

2. Image analysis in shear-SALS

The scattering pattern on the detector plate contains valuable information about the level of anisotropy and the average orientation angle of the glass fibers in the sample. This information is obtained using the following steps. First, the camera captures the scattering patterns as a function of time. Then, the background image is subtracted from each individual image to remove noise and the intensities originating solely from the primary laser beam. Next, the suspending rope of the beam stop is filtered from the image by interpolating the intensity level around the beam stop rope, resulting in a scattering pattern suitable for analysis, see Fig. 6(a). Then, intensity isocontours are computed (except for the beam stop region), which correspond to fixed q-values along the horizontal direction, through the beam stop (azimuthal coordinate of 90 °), which results in Fig. 6(b). Subsequently, ellipses are fitted through the points of the isocontour plot, see Fig. 6(c). The fitted ellipses are of the following form: A x 2 + B x y + C y 2 + D x + E y + F = 0. The eccentricity of these ellipses can be described using two eigenvalues, ν 1 = ( a / 2 ) 2 and ν 2 = ( b / 2 ) 2, where a and b are, respectively, the long and short axes of the ellipse. The anisotropy of the fibers according to the shear-SALS experiment, ξ SALS, is determined by the ratio of both eigenvalues, similar to Sec. II A,
(24)
The average angle of the glass fibers in the composite sample is equal to the orientation angle of the ellipse in the scattering pattern plus 90 ° and is determined using
(25)
where A, B, and C are the parameters describing the fitted ellipse. The angle of the ellipse corresponding to the long axis is considered so that ψ SALS is in the range [ π 4 , 3 π 4 ]. In the remainder of this work, the eccentricity will be determined from the ellipse fitted on the intermediate intensity values [red ellipse in Fig. 6(b)], as this ellipse effectively reduces interference from both the beam stop region and the detector edge. In Sec. S3 of the supplementary material, typical scattering patterns as a function of strain are shown, as well as the intensity as a function of q for different strains and for directions along the major and minor axes of the ellipse describing the intensity isocontours.
FIG. 6.

Image analysis procedure: (a) scattering pattern after background subtraction and after removal of the beamstop rope, (b) intensity isocontour plots, including ellipses describing the intensity isocontours, and (c) the eigenvalues and the orientation angle of the middle ellipse.

FIG. 6.

Image analysis procedure: (a) scattering pattern after background subtraction and after removal of the beamstop rope, (b) intensity isocontour plots, including ellipses describing the intensity isocontours, and (c) the eigenvalues and the orientation angle of the middle ellipse.

Close modal

In this section, the macroscopic shear field is validated by means of viscoelastic numerical simulations. First, the balance equations and boundary and initial conditions are discussed, afterward the numerical method is explained, and finally the macroscopic shear rate in the sample is determined. In these simulations, a cuboid of viscoelastic material is considered, where the velocity is prescribed on a part of the front and back surface, visualized with the gray area in Fig. 7. The gray area presents the part of the composite sample in the experimental setup that is attached to the clamps.

FIG. 7.

Schematic representation of the considered cuboid, including the important definitions to describe the geometry, boundary, and initial conditions. On the grayed surfaces, the velocity is prescribed. All remaining surfaces are denoted with Γ f. The curve on which the real shear rate is determined is shown as the orange line in the middle of the domain.

FIG. 7.

Schematic representation of the considered cuboid, including the important definitions to describe the geometry, boundary, and initial conditions. On the grayed surfaces, the velocity is prescribed. All remaining surfaces are denoted with Γ f. The curve on which the real shear rate is determined is shown as the orange line in the middle of the domain.

Close modal
The mass and momentum balance are used to describe the shear flow of the fluid. For an incompressible fluid, neglecting inertia and including gravity, the balance equations yield
(26)
(27)
where u is the fluid velocity vector, ρ is the polymer density ( ρ = 920 kg m 3), g = g e x the gravitation vector, with g being the gravitational constant ( g = 9.81 ms 2), and σ is the Cauchy stress tensor,
(28)
where p is the pressure, I is the unit tensor, η sol is the Newtonian solvent viscosity, and D = ( u + ( 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,
(29)
where η p is the polymer viscosity and λ is the polymer relaxation time. The latter two are related through the polymer modulus G = η p / λ. The conformation tensor is c, and its evolution equation for a Giesekus constitutive model is described by [44]
(30)
where α is a (dimensionless) constitutive parameter, which influences the shear thinning behavior and the symbol ( ) denotes the upper-convected time derivative,
(31)
where c ˙ is the material derivative.
The balance equations are closed using the boundary conditions and initial condition described in this section. The velocity on surfaces Γ 3 and Γ 4 is set to zero. On surfaces Γ 1 and Γ 2, the velocity in the x-direction is set to d γ ˙ a, with γ ˙ a being the applied shear rate, and the velocity in x- and z-directions is set to zero. All other surfaces, Γ f (see Fig. 7) are traction free. Summarizing, this gives
(32)
(33)
(34)
where n is the outwardly directed unit normal vector of surface Γ f.
Finally, a stress-free initial condition is assumed, i.e., at the initial time the viscoelastic stress is zero everywhere in the fluid,
(35)

The balance equations are solved using an in-house developed software package utilizing the finite element method (FEM). The elemental mesh is generated using Gmsh [45]. For the velocity and pressure, isoparametric tetrahedral P 2 P 1 (Taylor–Hood) elements are used, whereas for the velocity gradient and conformation, tetrahedral P 1 elements are used. Discrete Elastic-Viscous Stress Splitting using the G approach (DEVSS-G), Streamline Upwind Petrov Galerkin (SUPG), and the log representation of the conformation tensor are used to increase the stability of the simulations [46–50].

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 [51]. The second order semi-implicit Gear scheme was extended to not be limited in time step size by very fast modes, see our previous work [2]. After solving the conformation prediction, the nodal points in the mesh are updated using a Lagrangian Adams-Bashforth scheme. After a number of time steps, at large strains, the elements can become too distorted with possible inaccuracies as a result. However, the model prediction in Sec. II A shows that the orientation kinetics are fast and that most particles are already aligned in the velocity direction at low strains of γ 5. Therefore, a remeshing step is not implemented in these simulations. In Sec. S4 of the supplementary material, mesh and time step convergence tests are shown. According to these convergence tests, a strain step of γ ˙ a Δ t = 0.0055 is used with an element size equal to the cuboid thickness divided by 4, as this provides an optimum in computational efficiency and accuracy.

This section is divided into two main parts. First, characterization and verification results are shown. Here, the anisotropy of the in situ SALS images is compared with the anisotropy resulting from ex situ CT scans, and the local shear rates in the sample are analyzed by means of viscoelastic macroscopic simulations. In the second part, the in situ fiber orientation kinetics are investigated in shear flow. Afterward, the corresponding transient shear viscosity of the fiber composites is investigated. In this work, the fiber orientation kinetics and rheology are investigated on separate samples, with the rheology data covering a smaller strain range.

1. Comparison of in situ SALS with ex situ CT scans

In this section, the anisotropy of the fibers in the composites determined from the scattering images is compared to that resulting from the CT scans. Our previous study [1] showed via CT scans that most fibers in undeformed samples are oriented in the x y-plane. CT-scan samples are obtained as follows: First, shear-SALS experiments are performed using an applied shear rate of 0.03  s 1 until a strain of about 10 is reached, and the anisotropy ξ SALS is determined. Then, the setup, including the sample, is cooled with an external fan for about 10 min. Afterward, the sample is removed from the setup and cut into a size of 1 × 2 mm 2 at the location of the laser beam. Subsequently, a CT scan is performed to determine the anisotropy, ξ CT. The visualization in Geodict of a sheared CT-scanned sample having a fiber volume fraction of 1% is shown in Fig. 8. The visualizations in Geodict of sheared CT-scanned samples having fiber volume fractions of 5% and 10% are shown in Sec. S5 of the supplementary material. In Sec. S5 of the supplementary material, the visualization of CT-scanned samples sheared to a strain of about 3 is also shown for different fiber volume fractions. For all considered CT scans, the orientation component in the thickness direction is close to zero, meaning that the anisotropy can indeed be computed according to Eq. (21). To validate the orientation kinetics obtained using the SALS method, the anisotropy resulting from the scattering pattern is compared to that resulting from the CT scans. The ratio between the two is called the anisotropy ratio,
(36)
In general, the anisotropy ratio is dependent on the q-value and the fiber volume fraction. In this study, the q-range is slightly above 2 π / D. Since the q-value is fixed in this study, the anisotropy ratio is only computed for different volume fractions, see Table I, where both measurements are performed on three different samples to compute the average and the standard deviation. The table shows that for all considered fiber volume fractions, the anisotropy ratio is about 1, meaning that the SALS methodology is a reliable method to assess the anisotropy in the experimental system. In the remainder of this work, the corrected anisotropy, ξ = ξ SALS C ξ, will be considered, which results in an insignificant correction.
FIG. 8.

Representation in Geodict of a CT-scanned sample. The sample has a fiber volume fraction of 1% and is sheared with an applied rate of γ ˙ a = 0.03 s 1 until a strain of about 10.

FIG. 8.

Representation in Geodict of a CT-scanned sample. The sample has a fiber volume fraction of 1% and is sheared with an applied rate of γ ˙ a = 0.03 s 1 until a strain of about 10.

Close modal
TABLE I.

Ratio between anisotropy from x-ray tomography and anisotropy from the scattering images for different fiber volume fractions, after shear until an applied strain of 10.

1% Fibers5% Fibers10% Fibers
Cξ 0.99 ± 0.01 1.001 ± 0.001 1.01 ± 0.01 
ξCT 0.95 ± 0.005 0.97 ± 0.004 0.97 ± 0.008 
ξSALS 0.96 ± 0.01 0.97 ± 0.003 0.96 ± 0.01 
1% Fibers5% Fibers10% Fibers
Cξ 0.99 ± 0.01 1.001 ± 0.001 1.01 ± 0.01 
ξCT 0.95 ± 0.005 0.97 ± 0.004 0.97 ± 0.008 
ξSALS 0.96 ± 0.01 0.97 ± 0.003 0.96 ± 0.01 

2. Validation of shear rate in the sample at the laser location

To validate the local shear rates occurring in the composite samples in the shear-SALS setup, macroscopic viscoelastic shear simulations are performed, see Sec. IV. A seven-mode and a five-mode Giesekus material model is employed, describing the transient behavior of, respectively, LDPE and LLDPE at 150  °C. The Giesekus model parameters are determined in our previous work [2]. The solvent viscosity parameter, η sol = 0 Pa s, in these simulations, except in the first time step, where η sol = 1.0 Pa s for numerical reasons. The deformed shape of a sheared sample in the viscoelastic simulations is compared to the deformed shape of a real sample sheared in the shear-SALS setup in Fig. 9. Both samples in the figure are an LDPE sheared with an applied shear rate of γ a ˙ = 0.03 s 1. The difference between the left [Figs. 9(a) and 9(b)] and right [Figs. 9(c) and 9(d)] samples is the gap size, which is, respectively, 4.0 and 1.8 mm. The deformed shapes of the real and simulated samples agree well, which provides confidence in the simulations. In the simulations, the rate-of-deformation tensor is computed spatially as a function of time, which is required to compute the real shear rate ( γ ˙ r = 2 D x y). The real shear rate is computed in the middle of the x z cross section, ranging from the top of the sample (positive y-coordinate) to the bottom of the sample (negative y-coordinate). The distance from the top of the sample to the bottom of the sample will be referred to as the arc length and is shown as the orange curve in Fig. 10.

FIG. 9.

Comparison between [(a) and (c)] the shape of a sample after being sheared in the shear-SALS setup and [(b) and (d)] the deformed shape resulting from the macroscopic viscoelastic simulations. (a) and (b) show the comparison for a sample with an initial gap of 4.0 mm, and (c) and (d) show the comparison for a sample with an initial gap of 1.8 mm.

FIG. 9.

Comparison between [(a) and (c)] the shape of a sample after being sheared in the shear-SALS setup and [(b) and (d)] the deformed shape resulting from the macroscopic viscoelastic simulations. (a) and (b) show the comparison for a sample with an initial gap of 4.0 mm, and (c) and (d) show the comparison for a sample with an initial gap of 1.8 mm.

Close modal
FIG. 10.

Real shear rate, γ ˙ r, as a function of arc length (see Fig. 7) according to macroscopic shear simulations for an LDPE material model with a sample gap of 1.9 mm and a sample thickness of 1.5 mm at different shear strains. The horizontal dashed line indicates the applied shear rate, γ ˙ a, and the zone between the vertical dotted lines indicates the location of the fiber orientation measurements.

FIG. 10.

Real shear rate, γ ˙ r, as a function of arc length (see Fig. 7) according to macroscopic shear simulations for an LDPE material model with a sample gap of 1.9 mm and a sample thickness of 1.5 mm at different shear strains. The horizontal dashed line indicates the applied shear rate, γ ˙ a, and the zone between the vertical dotted lines indicates the location of the fiber orientation measurements.

Close modal
The real shear rate as a function of the arc length at different applied strains ( γ a = γ ˙ a t) is plotted in Fig. 10 using a seven-mode Giesekus material model representing the behavior of LDPE. In this viscoelastic simulation, the velocity on the moving surfaces is such that the applied shear rate is γ ˙ a = 0.55 s 1 (see horizontal dashed line in Fig. 10). The zone between the vertical dotted lines in the figure indicates the location where the laser beam passes through the composite sample and, thus, indicates the location of the fiber orientation measurement. At this position, the ratio between the real shear rate according to the viscoelastic simulations and the applied shear rate is defined as
(37)
Figure 10 shows that for strains γ a < 7, C γ ˙ 0.77, implying that the real strain in the composite sample at the laser location in the shear-SALS setup is lower than the applied strain. For strains γ a > 7, fibers are already oriented in the flow direction (see Sec. V B 1) and the exact value of C γ ˙ becomes irrelevant. Also, for strains γ a > 7, rheological measurements are not possible as explained in Sec. V B 2. Therefore, in the remainder of this work, C γ ˙ is taken as a constant per parameter set, which will be determined at low strains at the location of the fiber orientation measurement.

The ratio between the real shear rate and applied shear rate is a function of the sample thickness, sample gap, and the Weissenberg number and is shown for LDPE and LLDPE in Fig. 11. Each panel in Fig. 11 shows a contour plot of the shear rate ratio at a single applied shear rate, which is γ ˙ a = 0.03, 0.55, and 5 s 1 for the left, middle, and right panels, respectively. For LDPE [Figs. 11(a)11(c)], the shear rate ratio is strongly dependent on the applied rate because elastic stresses dominate the material behavior. For LLDPE [Figs. 11(d)11(f)], the influence of applied rate on the shear rate ratio is negligible, because for all applied rates there are limited elastic effects. In the remainder of this work, the real strain, γ r = t γ ˙ a C γ ˙, will be considered.

FIG. 11.

Dependence of the ratio between the real shear rate in the middle of the sample and the macroscopically applied shear rate on sample thickness, sample gap, and applied shear rate, according to viscoelastic simulations. (a), (b), and (c) correspond to LDPE, which is described using a seven-mode viscoelastic Giesekus material model. (d), (e), and (f) correspond to LLDPE, which is described using a five-mode viscoelastic Giesekus material model. In the left, middle, and right figures, the applied shear rates are γ ˙ a = 0.03, 0.55, and 5.0 s 1, respectively.

FIG. 11.

Dependence of the ratio between the real shear rate in the middle of the sample and the macroscopically applied shear rate on sample thickness, sample gap, and applied shear rate, according to viscoelastic simulations. (a), (b), and (c) correspond to LDPE, which is described using a seven-mode viscoelastic Giesekus material model. (d), (e), and (f) correspond to LLDPE, which is described using a five-mode viscoelastic Giesekus material model. In the left, middle, and right figures, the applied shear rates are γ ˙ a = 0.03, 0.55, and 5.0 s 1, respectively.

Close modal

1. Fiber orientation in shear

In this section, the fiber orientation kinetics are investigated for shear flow as a function of Weissenberg number, fiber volume fraction, and matrix rheology. The fiber orientation kinetics are measured using SALS as described in Sec. III D 2, and the results are shown in Figs. 12 and 13. In both figures, (a) and (c) show, respectively, the average anisotropy and the corresponding average orientation angle of LDPE, whereas (b) and (d) show these data for LLDPE. The solid black lines in the figures correspond to the average anisotropy and the corresponding average orientation angle as predicted by the multiparticle model, see Eqs. (15) and (16). The particle shape factor in the model corresponds to a prolate ellipsoid with an aspect ratio of 10, which is the same geometry as used in our study on orientation in uniaxial extension. On the x-axis, the real strain is plotted, which is the applied strain multiplied by the shear rate correction factor as discussed in Sec. V A 2. The solid markers in the graphs indicate averages with error bars indicating the standard deviation, which is computed using at least three repeat measurements. The initial anisotropy of the individual samples varied significantly, and the analysis only took samples into account with a minimum anisotropy of 0.5 or lower. The initial strain of these samples is corrected so that the minimum anisotropy and the model prediction overlap. Individual measurement data at strains lower than the strain corresponding to minimum anisotropy are plotted with open markers in the figures. Furthermore, the average orientation angle corresponding to low anisotropies is not well defined, which at small strains results in large standard deviations in the average orientation angle data.

FIG. 12.

[(a) and (b)] Average anisotropy and [(c) and (d)] the corresponding average orientation angle, for different Weissenberg numbers, as a function of shear strain for [(a) and (c)] LDPE and for [(b) and (d)] LLDPE as matrix material. All samples have a fiber volume fraction of 1%. For LDPE and LLDPE, the applied Weissenberg numbers range, respectively, from 1.6 to 265 and from 0.02 to 3.1. The black solid lines correspond to the prediction of the multiparticle model [Eqs. (15) and (16)]. The filled markers indicate averages with error bars indicating the standard deviation. The open markers indicate individual experimental data points.

FIG. 12.

[(a) and (b)] Average anisotropy and [(c) and (d)] the corresponding average orientation angle, for different Weissenberg numbers, as a function of shear strain for [(a) and (c)] LDPE and for [(b) and (d)] LLDPE as matrix material. All samples have a fiber volume fraction of 1%. For LDPE and LLDPE, the applied Weissenberg numbers range, respectively, from 1.6 to 265 and from 0.02 to 3.1. The black solid lines correspond to the prediction of the multiparticle model [Eqs. (15) and (16)]. The filled markers indicate averages with error bars indicating the standard deviation. The open markers indicate individual experimental data points.

Close modal
FIG. 13.

[(a) and (b)] Average anisotropy and [(c) and (d)] the corresponding average orientation angle, for different fiber volume fractions, as a function of shear strain for [(a) and (c)] LDPE and for [(b) and (d)] LLDPE as matrix material. All samples are sheared with an applied shear rate of γ ˙ a = 0.55 s 1. The black solid lines correspond to the prediction of the multiparticle model [Eqs. (15) and (16)]. The filled markers indicate averages with error bars indicating the standard deviation. The open markers indicate individual experimental data points.

FIG. 13.

[(a) and (b)] Average anisotropy and [(c) and (d)] the corresponding average orientation angle, for different fiber volume fractions, as a function of shear strain for [(a) and (c)] LDPE and for [(b) and (d)] LLDPE as matrix material. All samples are sheared with an applied shear rate of γ ˙ a = 0.55 s 1. The black solid lines correspond to the prediction of the multiparticle model [Eqs. (15) and (16)]. The filled markers indicate averages with error bars indicating the standard deviation. The open markers indicate individual experimental data points.

Close modal

The influence of Weissenberg number on the fiber orientation kinetics for composites with a fiber volume fraction of 1% is shown in Fig. 12. For both LDPE and LLDPE, the applied shear rate in the experiments ranged from γ ˙ = 0.03 to 5.0  s 1, which corresponds to an applied viscosity-averaged Weissenberg number ranging from Wi = 1.6–265 for LDPE and Wi = 0.02–3.1 for LLDPE. The initial anisotropies for the samples sheared at the highest shear rate are higher because of the negative prestrain, as explained inSec. III D. From Fig. 12, the first key result of this work can be obtained: at strains below 10, the fiber orientation kinetics are independent of the Weissenberg number. For both materials with Weissenberg numbers ranging over several decades, the confidence interval of the experimental measurements overlaps with the model prediction. Upon closer examination, it is observed that the average anisotropy of LDPE composites subjected to the two highest shear rates is slightly below the model prediction at a strain of 3–5, indicating a slightly slower orientation.

The influence of fiber volume fraction on the fiber orientation kinetics is investigated in Fig. 13. All composite samples are sheared with an applied shear rate γ ˙ a = 0.55 s 1. The figure shows the second key result of this work. For concentrations in the dilute, semidilute, and up to the onset of the concentrated regime ( ϕ D / L), there is no influence of fiber volume fraction on the fiber orientation kinetics, for these short glass fibers, since the confidence intervals of the different measurements overlap. Furthermore, it is observed that the experimentally obtained fiber orientation kinetics are reasonably well described by the model. The reader is reminded that the particle kinetics in the model are based on Jeffery’s equation and do not include any orientation kinetics reduction parameter (e.g., strain reduction factor, Brownian kinetics parameter, that are used in [9,10,21,23–25]). This remarkable result will be further discussed in Sec. VI. In Sec. S6 of the supplementary material, the fiber orientation kinetics are shown for LDPE composites at different fiber volume fractions and different Weissenberg numbers. The results confirm that both variables do not influence the fiber orientation kinetics.

2. Shear rheology of fiber composites

The transient shear rheology of the fiber composites is measured according to the method described in Sec. III D 1 using the same shear-SALS setup and compared to the rheological model prediction of Eq. (20). The samples considered in this section have a thickness varying between 1 and 1.4 mm, and the fibers are distributed over the complete thickness. This is in contrast to the 5% and 10% composite samples in Sec. V B 1, where the fiber-filled thickness of the samples was, respectively, 0.44 and 0.22 mm. The transient shear viscosity of fiber composites with initially an isotropic fiber orientation distribution as a function of time is shown for LDPE and LLDPE in, respectively, Figs. 14(a) and 14(b). The transient shear viscosity is also measured for LDPE and LLDPE composites presheared to a strain of γ r , pre 2.0, see Figs. 15(a) and 15(b), respectively. The x-axis on the top of the plot indicates the considered strains for the data from the sliding plate measurements. Before conducting the transient shear viscosity measurements of the presheared samples, relaxation of the stress in the material was ensured. According to Figs. 12 and 13, the presheared composite samples are expected to have an initial anisotropy of ξ 0.7 and an initial average fiber orientation angle of ψ 65 °. In Figs. 14 and 15, the estimated values for the anisotropy are included in the graph. A steady state fiber orientation is not obtained in these experiments. The presheared composites in Fig. 15 are expected to have an isotropic orientation after about 4 s.

FIG. 14.

Transient shear viscosity of the fiber composites as a function of time. (a) LDPE sheared with γ ˙ a = 0.55 s 1 ( γ ˙ r 0.43 s 1) and (b) LLDPE sheared with γ ˙ a = 0.55 s 1 ( γ ˙ r 0.43 s 1) for samples with initially an isotropic fiber orientation ( ξ 0) evolving to semianisotropic ( ξ 0.8, ψ 70 °). The x-axis on the top of the plot indicates the considered strains for the data from the sliding plate measurements. The error bars indicate the standard deviation. The dashed line corresponds to the transient shear viscosity of the unfilled matrix, measured in standard rotational mode. The solid lines correspond to the rheological model prediction [Eq. (20)]. The white area indicates the region where accurate data are obtained.

FIG. 14.

Transient shear viscosity of the fiber composites as a function of time. (a) LDPE sheared with γ ˙ a = 0.55 s 1 ( γ ˙ r 0.43 s 1) and (b) LLDPE sheared with γ ˙ a = 0.55 s 1 ( γ ˙ r 0.43 s 1) for samples with initially an isotropic fiber orientation ( ξ 0) evolving to semianisotropic ( ξ 0.8, ψ 70 °). The x-axis on the top of the plot indicates the considered strains for the data from the sliding plate measurements. The error bars indicate the standard deviation. The dashed line corresponds to the transient shear viscosity of the unfilled matrix, measured in standard rotational mode. The solid lines correspond to the rheological model prediction [Eq. (20)]. The white area indicates the region where accurate data are obtained.

Close modal
FIG. 15.

Transient shear viscosity of the fiber composites as a function of time. (a) LDPE sheared with γ ˙ a = 0.55 s 1 ( γ ˙ r 0.43 s 1) and (b) LLDPE sheared with γ ˙ a = 0.55 s 1 ( γ ˙ r 0.43 s 1) for presheared samples with initially a semianisotropic fiber orientation ( ξ 0.7, ψ 65 °) evolving to isotropic ( ξ 0) and back to semianisotropic ( ξ 0.7, ψ 65 °). The x-axis on the top of the plot indicates the considered strains for the data from the sliding plate measurements. The error bars indicate the standard deviation. The dashed line corresponds to the transient shear viscosity of the unfilled matrix, measured in standard rotational mode. The solid lines correspond to the rheological model prediction [Eq. (20)]. The white area indicates the region where accurate data are obtained.

FIG. 15.

Transient shear viscosity of the fiber composites as a function of time. (a) LDPE sheared with γ ˙ a = 0.55 s 1 ( γ ˙ r 0.43 s 1) and (b) LLDPE sheared with γ ˙ a = 0.55 s 1 ( γ ˙ r 0.43 s 1) for presheared samples with initially a semianisotropic fiber orientation ( ξ 0.7, ψ 65 °) evolving to isotropic ( ξ 0) and back to semianisotropic ( ξ 0.7, ψ 65 °). The x-axis on the top of the plot indicates the considered strains for the data from the sliding plate measurements. The error bars indicate the standard deviation. The dashed line corresponds to the transient shear viscosity of the unfilled matrix, measured in standard rotational mode. The solid lines correspond to the rheological model prediction [Eq. (20)]. The white area indicates the region where accurate data are obtained.

Close modal

All composite samples in Figs. 14 and 15 are sheared with an applied shear rate, γ a ˙ = 0.55 s 1, and in these figures, the experimentally obtained average transient shear viscosities are indicated with markers and the corresponding standard deviations with error bars, which are determined using at least three repeat measurements. The solid lines in the figures correspond to the transient shear viscosity predicted by the rheological fiber-composite model [Eq. (20)] which assumes a constant, random orientation of fibers. The dashed line corresponds to the transient shear viscosity of unfilled LDPE or LLDPE, measured using a conventional ARES rotational rheometer, see our previous work [1]. At small times, the transient shear viscosity of unfilled LDPE and LLDPE measured using the shear-SALS setup highly agrees with that measured using the ARES rotational rheometer. At larger times, the axial force signal drops significantly, because of the large free sample surface and the viscosity measurement becomes inaccurate, which is shown in the gray area in the figures. The drop of the axial force signal occurs at larger times for the presheared samples (Fig. 15) compared to the non-presheared samples (Fig. 14). At the start of the gray area, the composite samples are, according to Fig. 13, expected to have an anisotropy of ξ 0.8 and an average orientation angle of ψ 70 °. Due to this limitation in strain, this setup is only capable of measuring the steady state viscosity of viscous materials. For materials with longer relaxation times, the strain that can be achieved is insufficient for an accurate measurement of the steady state viscosity.

From Figs. 14 and 15, one key result can be obtained. The fiber composite rheological model [Eq. (20)], which assumes a random and isotropic fiber orientation, agrees well with the transient shear viscosity of the composite samples at different strains, thus different anisotropies ( ξ 0 to ξ 0.8 with ψ 70 °). Note, composites with 10% fiber volume fraction show a subtle local minimum when the fibers have an isotropic orientation distribution, see Figs. 15(a) and 15(b). In these figures, initially and at times just before the gray area, the semianisotropy and corresponding average orientation angle of the fibers is, respectively, aligned in the compression and the extension direction of the shear flow (see the text in the figures), increasing the composites’ viscosity [11]. The influence of solely the fibers on the composites’ viscosity is shown in more detail by rescaling the data with the matrix viscosity, see Sec. S7 of the supplementary material.

The initial fiber orientation kinetics after startup of shear flow for suspensions with concentrations up to the concentrated regime, ϕ D / L, presented in this work (Fig. 13) can be described using Jeffery’s equation. This observation is contradictory to most of the existing literature, where the kinetics are described using the extended Folgar–Tucker model, which includes phenomenological terms to slow down the orientation kinetics. One such term is the Brownian diffusion term ( C I), and another term is the slip factor or strain reduction factor ( κ) [11]. In Table II, a comparison is made between the phenomenological model parameters to describe the fiber orientation kinetics in shear flow presented in this work and that of the existing literature. The table also includes details about the corresponding experimental systems. The “slip factor” in the majority of the existing literature ranges between 0.25 and 0.33, implying that the orientation kinetics are three to four times slower than predicted using Jeffery’s equation. A possible explanation is that the fiber aspect ratio in these systems is larger, generating higher stresses on the fibers, which increases the possibility that fibers are bending and this might affect the orientation kinetics. Fiber buckling in shear flow can be estimated by the equation of Forgacs and Mason [54] and is not expected in the experiments considered in the current study.

TABLE II.

Experimental systems and corresponding Folgar–Tucker parameters to describe the orientation kinetics in pure shear flow. The κ and CI are the phenomenological “slip factor” or “strain reduction factor” and, respectively, the phenomenological Brownian diffusion factor used in Folgar–Tucker models to describe the orientation kinetics [11].

Sourceϕ (vol. %)reϕ ⋅ reD (μm)Matrix materialκCI
Eberle et al. [2115 28 4.2 12.9 Polybutylene terephthalate 0.3 0.006 
Cieslinski et al. [2314 52 7.28 13.3 Polypropylene 0.25 0.005 
Ortman et al. [1014 202 28.28 14.5 Polypropylene 0.25 0.005 
Ortman et al. [10202 10.10 14.5 Polypropylene 0.32 0.0035 
Sepehr et al. [91.5 20 0.2 14 Polybutene 0.31 0.002 
Sepehr et al. [93.2 20 0.64 14 Polybutene 0.26 0.0045 
Sepehr et al. [97.06 20 1.41 14 Polybutene 0.24 0.0055 
Sepehr et al. [2514 18.6 2.60 14 Polypropylene 0.33 0.0035 
Sepehr et al. [2510 18.6 1.77 14 Polypropylene 0.38 0.0067 
Wang et al. [24,5215 25 3.75 16.7 Polybutylene terephthalate 0.11 0.0112 
Férec et al. [533.8 20 0.76 16 Polybutene 0.995b 0.0009a,b 
Férec et al. [535.8 20 1.16 16 Polybutene 0.992b 0.0014a,b 
Férec et al. [538.1 20 1.62 16 Polybutene 0.988b 0.002a,b 
Férec et al. [5310.5 20 2.1 16 Polybutene 0.985b 0.0026a,b 
Current work 10 11 1.1 10 Polyethylene 
Sourceϕ (vol. %)reϕ ⋅ reD (μm)Matrix materialκCI
Eberle et al. [2115 28 4.2 12.9 Polybutylene terephthalate 0.3 0.006 
Cieslinski et al. [2314 52 7.28 13.3 Polypropylene 0.25 0.005 
Ortman et al. [1014 202 28.28 14.5 Polypropylene 0.25 0.005 
Ortman et al. [10202 10.10 14.5 Polypropylene 0.32 0.0035 
Sepehr et al. [91.5 20 0.2 14 Polybutene 0.31 0.002 
Sepehr et al. [93.2 20 0.64 14 Polybutene 0.26 0.0045 
Sepehr et al. [97.06 20 1.41 14 Polybutene 0.24 0.0055 
Sepehr et al. [2514 18.6 2.60 14 Polypropylene 0.33 0.0035 
Sepehr et al. [2510 18.6 1.77 14 Polypropylene 0.38 0.0067 
Wang et al. [24,5215 25 3.75 16.7 Polybutylene terephthalate 0.11 0.0112 
Férec et al. [533.8 20 0.76 16 Polybutene 0.995b 0.0009a,b 
Férec et al. [535.8 20 1.16 16 Polybutene 0.992b 0.0014a,b 
Férec et al. [538.1 20 1.62 16 Polybutene 0.988b 0.002a,b 
Férec et al. [5310.5 20 2.1 16 Polybutene 0.985b 0.0026a,b 
Current work 10 11 1.1 10 Polyethylene 
a

Maximum value, CI decreases as a function of orientation.

b

Conversion of parameters used by Férec et al. [53], see Sec. S8 in the supplementary material.

The majority of the systems considered in the existing literature are further in the concentrated regime, which can be seen by scaling the volume fraction with the inverse of the fiber aspect ratio, ϕ / ( D / L ) = ϕ r e (note that for the concentrated regime ϕ > D / L, or ϕ r e > 1). Exceptions are the work of Sepehr et al. [9] and Férec et al. [53], both of them determined the fiber orientation kinetics using an indirect plate-plate rheology method. First, they applied a negative preshear to promote tangential fiber alignment. After shearing in the opposite direction, at a certain strain, the viscosity is increased, which is then related to the orientation kinetics (tumbling) of the fibers. They both, thus, do not consider the orientation kinetics of samples with an initially random fiber orientation. To describe the orientation kinetics for composites in the semiconcentrated regime, ( D / L ) 2 ϕ D / L, Sepehr et al. [9] required a κ = 0.31 in the Folgar–Tucker model and Férec et al. [53] included a fiber-fiber interaction tensor in the fiber orientation evolution, a κ 1 (see Sec. S8 of the supplementary material for the conversion of the model parameters of their model to that of the extended Folgar–Tucker model). The parameters found by the former authors for an initially anisotropic orientation do not agree with the parameters presented in the current paper for an initially isotropic orientation, while the parameters found by the latter authors do. For Sepehr et al. [9] and Férec et al. [53], small differences in average orientation angle after preshear can significantly influence the strain at which the particles tumble (see Fig. 2) and, thus, the corresponding model parameters. Unfortunately, the exact average orientation angle after the preshear is unknown. Moreover, our results suggest that the fiber orientation kinetics are solely influenced by the applied strain. In the parallel plate geometry, the applied strain is a function of the radial coordinate, which might lead to uneven tumbling of the fibers with additional fiber-fiber interactions. This might result in inaccurate rheological measurements and hence an inaccurate determination of the fiber orientation kinetics. Nevertheless, the torque of the rotational rheometer scales with the cube of the radial coordinate, meaning that the viscosity is primarily influenced by the particles at larger radial positions, which undergo the applied shear rate and strain. Finally, in the work of Bounoua et al. [34] and Li et al. [35], the slower orientation kinetics were partially ascribed to confinement effects occurring in between parallel plates in the rotational rheometer. The confinement effects can also introduce shear banding phenomena, as observed by Li et al. [36], which also affects the fiber orientation kinetics.

The slower orientation kinetics of the experimental systems in the existing literature might also be explained by the occurrence of slip between the fiber surface and the matrix. To investigate the occurrence of slip in the experimental system considered in the current work, fiber orientation kinetics experiments are also performed for composites consisting of amino-silane-coated glass fibers in an LDPE-g-MAH matrix. Slip in these systems is prevented because the amino-silane coating can covalently bond to the maleic anhydride group of the matrix. Fiber orientation experiments for composites with LDPE-g-MAH matrices are performed at 140  °C because the transient rheological properties of LDPE-g-MAH at this temperature are similar to those of LDPE at 150  °C (see Sec. S9 of the supplementary material). This prevents the effects of the material rheology and isolates the effect of slip on the fiber orientation kinetics. SEM images were taken (see Sec. S10 of the supplementary material), in order to verify the connection between the amino-silane coating of the fiber surface and the maleic anhydride group of the LDPE-g-MAH. These images are, however, inconclusive. A comparison between the SEM images with an LDPE or LDPE-g-MAH matrix shows that it is likely that both matrices adhere well to the fiber surface. In Fig. 16, a comparison is shown for the fiber orientation kinetics for glass fiber composites with an LDPE and an LDPE-g-MAH matrix. The considered composites both contained 1 vol. % glass fibers and were sheared with a rate of γ ˙ a = 0.55 s 1. Figures 16(a) and 16(b) represent the average anisotropy and average orientation angle, respectively. The x-axis, filled markers, open markers, and solid lines represent the same as explained in Figs. 12 and 13. The average anisotropy and average orientation angle for both composites overlap, implying that either slip did not occur in the composites without maleic anhydride groups or that the presence of slip does not influence the fiber orientation kinetics at these small strains.

FIG. 16.

(a) Anisotropy and (b) average orientation angle as a function of shear strain for LDPE and LDPE-g-MAH with 1% fiber volume fraction. In all experiments, the applied shear rate, γ ˙ a = 0.55 s 1. The solid black lines correspond to the prediction of the multiparticle model [Eqs. (15) and (16)]. The filled markers indicate averages with error bars indicating the standard deviation. The open markers indicate individual experimental data points.

FIG. 16.

(a) Anisotropy and (b) average orientation angle as a function of shear strain for LDPE and LDPE-g-MAH with 1% fiber volume fraction. In all experiments, the applied shear rate, γ ˙ a = 0.55 s 1. The solid black lines correspond to the prediction of the multiparticle model [Eqs. (15) and (16)]. The filled markers indicate averages with error bars indicating the standard deviation. The open markers indicate individual experimental data points.

Close modal

Finally, in the present study, the considered strains are limited and as a result, the stress in the material is not fully developed. Steady state stresses are, therefore, not reached. The fiber orientation in a developed stress field can still be dependent on the Weissenberg number or differ from Jeffery’s equation. Though, Fig. 12 also shows, indirectly, that a certain level of stress in the matrix does not influence the fiber orientation kinetics. The experiments performed with an applied shear rate of γ ˙ a = 5 s 1 are presheared to a negative strain. After applying shear in the positive direction, the minimum anisotropy is obtained after the prestrain amount and here already stress is present in the matrix. The resulting orientation kinetics overlap with the kinetics obtained for the other shear rates ( γ ˙ a = 0.55 s 1 and γ ˙ a = 0.03 s 1), where minimum anisotropy is obtained at a negligible stress level in the matrix.

In this study, the orientation kinetics and the corresponding rheology of molten short fiber-polymer composites were experimentally investigated for shear flow. A commercial rotational rheometer was modified with custom components to apply shear flow in the sliding-plate configuration and to measure the fiber orientation kinetics in the velocity–velocity gradient plane in situ using SALS. Moreover, the transient shear viscosity is determined using the standard axial force data of the commercial rheometer. The shear field in the composite during the experiments was validated through viscoelastic simulations, which indicate that the applied shear rate must be corrected with a constant factor to accurately determine the real shear rate in the sample at the location of the fiber orientation measurement. This correction factor is dependent on the sample thickness, the sample gap, and the Weissenberg number. The initial fiber orientation kinetics after startup of shear flow for non-Brownian short glass fibers suspended in different polyethylenes were studied as a function of Weissenberg number and fiber volume fraction. The Weissenberg number was varied over three decades, by varying the shear rate, and the fiber orientation kinetics remained unaffected. Furthermore, for concentrations in the dilute, semidilute, and up to the onset of the concentrated regime ( ϕ = D / L) no influence of fiber volume fraction on the fiber orientation kinetics was found for these short glass fibers. These orientation kinetics are reasonably well described by a model based on Jeffery’s equation, without the need to introduce additional terms to account for fiber-fiber interactions. The period between tumbling was not investigated in this study. The experimental setup established in the current work can, in the future, be used to investigate the orientation kinetics of other systems like long glass fiber composites. Moreover, in ongoing work, the experimental results are used to validate finite element simulations of fiber composites, which can then be used to predict fiber orientation in processing flows.

Finally, the transient shear viscosity of short fiber-filled composites was investigated for small strains ( γ < 2.6), where the fiber orientation remains isotropic to semianisotropic. For the considered system and these small strains, the viscosity is well described by a rheological model, which also assumes a constant isotropic fiber orientation.

See supplementary material for the comparison between the multiparticle model with a fully random initial orientation of particles and with an initial random particle orientation in the velocity–velocity gradient plane. Here also the calculation of the Rouse and reptation based Weisenberg numbers is shown. Moreover, typical scattering patterns as a function of shear strain are shown. Additionally, mesh and time convergence tests for the simulation of Sec. IV are shown. Moreover, ex situ x-ray tomography scans of sheared fiber composites are visualized. This includes additional fiber orientation kinetics measurements. Also, the influence of the fibers on the composites’ viscosity is shown. It also shows a conversion for model parameters between two different fiber orientation models. Moreover, a comparison between the transient rheology of LDPE and that of LDPE-g-MAH is shown. Finally, the supplementary material includes scanning electron microscopy images showing fiber-matrix adhesion for LDPE and LDPE-g-MAH composites.

This research forms part of the research programme of Dutch Polymer Institute, project #840 Engineering the rheology ANd processinG-induced structural anisotropy of poLymEr composites with non-Brownian fibrous particles (ANGLE). The authors thank Nippon Electric Glass NL for supplying the short glass fibers and Frank Huijnen (The Compound Company) for supplying the LDPE-g-MAH. Furthermore, the authors thank Hamid Ahmadi for preparing the samples for SEM and Marc van Maris for performing the SEM imaging. Finally, the authors thank M. A. Hulsen at the Eindhoven University of Technology (TU/e) for access to the TFEM software libraries.

The authors have no conflicts to disclose.

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

1.
Egelmeers
,
T. R. N.
,
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
(
2
),
171
185
(
2024
).
2.
Egelmeers
,
T. R. N.
,
R. M.
Cardinaels
,
P. D.
Anderson
, and
N. O.
Jaensson
, “
Numerical simulation of fiber orientation kinetics and rheology of fiber-filled polymers in uniaxial extension
,”
Phys. Fluids
36
(
2
),
023356
(
2024
).
3.
Jeffery
,
G. B.
, “
The motion of ellipsoidal particles immersed in a viscous fluid
,”
Proc. R. Soc. A
102
,
161
179
(
1922
).
4.
Folgar
,
F.
, and
C. L.
Tucker
, “
Orientation behavior of fibers in concentrated suspensions
,”
J. Reinf. Plast. Compos.
3
,
98
119
(
1984
).
5.
Eberle
,
A. P. R.
,
D. G.
Baird
, and
P.
Wapperom
, “
Rheology of non-Newtonian fluids containing glass fibers: A review of experimental literature
,”
Ind. Eng. Chem. Res.
47
(
10
),
3470
3488
(
2008
).
6.
Petrich
,
M. P.
,
D. L.
Koch
, and
C.
Cohen
, “
An experimental determination of the stress–microstructure relationship in semi-concentrated fiber suspensions
,”
J. Non-Newtonian Fluid Mech.
95
(
2-3
),
101
133
(
2000
).
7.
Pujari
,
S.
,
S.
Rahatekar
,
J. W.
Gilman
,
K. K.
Koziol
,
A. H.
Windle
, and
W. R.
Burghardt
, “
Shear-induced anisotropy of concentrated multiwalled carbon nanotube suspensions using x-ray scattering
,”
J. Rheol.
55
(
5
),
1033
1058
(
2011
).
8.
Sepehr
,
M.
,
P. J.
Carreau
,
M.
Moan
, and
G.
Ausias
, “
Rheological properties of short fiber model suspensions
,”
J. Rheol.
48
(
5
),
1023
1048
(
2004
).
9.
Sepehr
,
M.
,
P. J.
Carreau
,
M.
Grmela
,
G.
Ausias
, and
P. G.
Lafleur
, “
Comparison of rheological properties of fiber suspensions with model predictions
,”
J. Polym. Eng.
24
(
6
),
579
610
(
2004
).
10.
Ortman
,
K.
,
D. G.
Baird
,
P.
Wapperom
, and
A.
Whittington
, “
Using startup of steady shear flow in a sliding plate rheometer to determine material parameters for the purpose of predicting long fiber orientation
,”
J. Rheol.
56
(
4
),
955
981
(
2012
).
11.
Tucker
,
C. L.
,
Fundamentals of Fiber Orientation Description, Measurement and Prediction
(
Carl Hanser Verlag
,
Munchen
,
2023
).
12.
Deboeuf
,
S.
,
N.
Lenoir
,
D.
Hautemayou
,
M.
Bornert
,
F.
Blanc
, and
G.
Ovarlez
, “
Imaging non-Brownian particle suspensions with x-ray tomography: Application to the microstructure of Newtonian and viscoplastic suspensions
,”
J. Rheol.
62
(
2
),
643
663
(
2018
).
13.
Pujari
,
S.
,
L.
Dougherty
,
C.
Mobuchon
,
P. J.
Carreau
,
M. C.
Heuzey
, and
W. R.
Burghardt
, “
X-ray scattering measurements of particle orientation in a sheared polymer/clay dispersion
,”
Rheol. Acta
50
,
3
16
(
2011
).
14.
Dykes
,
L. M. C.
,
J. M.
Torkelson
, and
W. R.
Burghardt
, “
Shear-induced orientation in well-exfoliated polystyrene/clay nanocomposites
,”
Macromolecules
45
(
3
),
1622
1630
(
2012
).
15.
Fry
,
D.
,
B.
Langhorst
,
H.
Wang
,
M. L.
Becker
,
B. J.
Bauer
,
E. A.
Grulke
, and
E. K.
Hobbie
, “
Rheo-optical studies of carbon nanotube suspensions
,”
J. Chem. Phys.
124
(
5
),
054703
(
2006
).
16.
Mehdipour
,
M. N.
,
N.
Reddy
,
R. J.
Shor
, and
G.
Natale
, “
Orientation dynamics of anisotropic and polydisperse colloidal suspensions
,”
Phys. Fluids
34
(
8
),
083317
(
2022
).
17.
Lettinga
,
M. P.
, and
J. K. G.
Dhont
, “
Non-equilibrium phase behaviour of rod-like viruses under shear flow
,”
J. Phys.: Condens. Matter
16
(
38
),
S3929
S3939
(
2004
).
18.
Lettinga
,
M. P.
,
Z.
Dogic
,
H.
Wang
, and
J.
Vermant
, “
Flow behavior of colloidal rodlike viruses in the nematic phase
,”
Langmuir
21
(
17
),
8048
8057
(
2005
).
19.
Kang
,
K.
,
M. P.
Lettinga
,
Z.
Dogic
, and
J. K. G.
Dhont
, “
Vorticity banding in rodlike virus suspensions
,”
Phys. Rev. E: Stat. Nonlinear Soft Matter Phys.
74
(
2
),
026307
(
2006
).
20.
Stover
,
C. A.
,
D. L.
Koch
, and
C.
Cohen
, “
Observations of fibre orientation in simple shear flow of semi-dilute suspensions
,”
J. Fluid Mech.
238
,
277
296
(
1992
).
21.
Eberle
,
A. P. R.
,
G. M.
Vélez-García
,
D. G.
Baird
, and
P.
Wapperom
, “
Fiber orientation kinetics of a concentrated short glass fiber suspension in startup of simple shear flow
,”
J. Non-Newtonian Fluid Mech.
165
(
3-4
),
110
119
(
2010
).
22.
Cieslinski
,
M. J.
,
D. G.
Baird
, and
P.
Wapperom
, “
Obtaining repeatable initial fiber orientation for the transient rheology of fiber suspensions in simple shear flow
,”
J. Rheol.
60
(
1
),
161
174
(
2016
).
23.
Cieslinski
,
M. J.
,
P.
Wapperom
, and
D. G.
Baird
, “
Fiber orientation evolution in simple shear flow from a repeatable initial fiber orientation
,”
J. Non-Newtonian Fluid Mech.
237
,
65
75
(
2016
).
24.
Wang
,
J.
,
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
).
25.
Sepehr
,
M.
,
G.
Ausias
, and
P. J.
Carreau
, “
Rheological properties of short fiber filled polypropylene in transient shear flow
,”
J. Non-Newtonian Fluid Mech.
123
(
1
),
19
32
(
2004
).
26.
Laun
,
H. M.
, “
Orientation effects and rheology of short glass fiber-reinforced thermoplastics
,”
Colloid Polym. Sci.
262
,
257
269
(
1984
).
27.
Ganani
,
E.
, and
R. L.
Powell
, “
Rheological properties of rodlike particles in a Newtonian and a non-Newtonian fluid
,”
J. Rheol.
30
(
5
),
995
1013
(
1986
).
28.
Ramazani
,
S. A.
,
A.
Ait-Kadi
, and
M.
Grmela
, “
Rheology of fiber suspensions in viscoelastic media: Experiments and model predictions
,”
J. Rheol.
45
(
4
),
945
962
(
2001
).
29.
Ivanov
,
Y. A.
,
T. G. M.
Van de Ven
, and
S. G.
Mason
, “
Damped oscillations in the viscosity of suspensions of rigid rods. I. Monomodal suspensions
,”
J. Rheol.
26
(
2
),
213
230
(
1982
).
30.
Marín-Santibáñez
,
B. M.
,
J.
Pérez-González
, and
L.
de Vargas
, “
Shear rheometry and visualization of glass fiber suspensions
,”
Rheol. Acta
49
,
177
189
(
2010
).
31.
Ortman
,
K.
,
N.
Agarwal
,
A. P. R.
Eberle
,
D. G.
Baird
,
P.
Wapperom
, and
A. J.
Giacomin
, “
Transient shear flow behavior of concentrated long glass fiber suspensions in a sliding plate rheometer
,”
J. Non-Newtonian Fluid Mech.
166
(
16
),
884
895
(
2011
).
32.
Ericsson
,
K. A.
,
S.
Toll
, and
J.-A. E.
Månson
, “
Sliding plate rheometry of planar oriented concentrated fiber suspension
,”
Rheol. Acta
36
,
397
405
(
1997
).
33.
Laun
,
H. M.
, “
Elastic properties of polyethylene melts at high shear rates with respect to extrusion
,”
Rheol. Acta
21
,
464
469
(
1982
).
34.
Bounoua
,
S. N.
,
P.
Kuzhir
, and
E.
Lemaire
, “
Shear reversal experiments on concentrated rigid fiber suspensions
,”
J. Rheol.
63
(
5
),
785
798
(
2019
).
35.
Li
,
B.
,
Y.
Guo
,
P.
Steeman
,
M.
Bulters
, and
W.
Yu
, “
Wall effect on the rheology of short-fiber suspensions under shear
,”
J. Rheol.
65
(
6
),
1169
1185
(
2021
).
36.
Li
,
B.
,
W.
You
,
S.
Liu
,
L.
Peng
,
X.
Huang
, and
W.
Yu
, “
Role of confinement in the shear banding and shear jamming in noncolloidal fiber suspensions
,”
Soft Matter
19
(
46
),
8965
8977
(
2023
).
37.
Hinch
,
E. J.
, 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
).
38.
Dinh
,
S. M.
, and
R. C.
Armstrong
, “
A rheological equation of state for semiconcentrated fiber suspensions
,”
J. Rheol.
28
(
3
),
207
227
(
1984
).
39.
Hinch
,
E. J.
, and
L. G.
Leal
, “
The effect of Brownian motion on the rheological properties of a suspension of non-spherical particles
,”
J. Fluid Mech.
52
(
4
),
683
712
(
1972
).
40.
Lipscomb
,
G. G.
,
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
).
41.
Giesekus
,
H.
, “
Elasto-viskose flüssigkeiten, für die in stationären schichtströmungen sämtliche normalspannungskomponenten verschieden groß sind
,”
Rheol. Acta
2
,
50
62
(
1962
).
42.
Tucker
,
C. L.
, “
Flow regimes for fiber suspensions in narrow gaps
,”
J. Non-Newtonian Fluid Mech.
39
(
3
),
239
268
(
1991
).
43.
Math2Market GmbH, Germany, Geodict simulation software release 2022.
44.
Giesekus
,
H.
, “
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
).
45.
Geuzaine
,
G.
, and
J. F.
Remacle
, “
Gmsh: A 3D finite element mesh generator with built-in pre- and post-processing facilities
,”
Int. J. Numer. Methods Eng.
79
(
11
),
1309
1331
(
2009
).
46.
Guénette
,
R.
, and
M.
Fortin
, “
A new mixed finite element method for computing viscoelastic flows
,”
J. Non-Newtonian Fluid Mech.
60
(
1
),
27
52
(
1995
).
47.
Bogaerds
,
A. C. B.
,
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
).
48.
Brooks
,
A. N.
, 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
).
49.
Hulsen
,
M. A.
,
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
).
50.
Fattal
,
R.
, and
R.
Kupferman
, “
Constitutive laws for the matrix-logarithm of the conformation tensor
,”
J. Non-Newtonian Fluid Mech.
123
(
2
),
281
285
(
2004
).
51.
D’Avino
,
G.
, 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
).
52.
Chu
,
J.
, and
J. L.
Sullivan
, “
Recyclability of a glass fiber poly(butylene terephthalate) composite
,”
Polym. Compos.
17
(
3
),
523
531
(
1996
).
53.
Férec
,
J.
,
G.
Ausias
,
M. C.
Heuzey
, and
P. J.
Carreau
, “
Modeling fiber interactions in semiconcentrated fiber suspensions
,”
J. Rheol.
53
(
1
),
49
72
(
2009
).
54.
Forgacs
,
O. L.
, and
S. G.
Mason
, “
Particle motions in sheared suspensions: IX. Spin and deformation of threadlike particles
,”
J. Colloid Sci.
14
(
5
),
457
472
(
1959
).