1. Introduction
Particle-filled polymeric or viscoelastic liquids are widely used in industrial applications such as hydraulic fracturing (Barbati et al. Reference Barbati, Desroches, Robisson and McKinley2016), fibre spinning and extrusion moulding (Nakajima et al. Reference Nakajima, Kajiwara and McIntyre1994; Breitenbach Reference Breitenbach2002; Huang et al. Reference Huang, Zhang, Kotaki and Ramakrishna2003). In hydraulic fracturing, polymeric liquids are chosen for their high viscosity, and spherical particles are added as proppants to keep fractures open. In extrusion moulding and fibre spinning, solid particles enhance the strength of the final product. These processes often involve uniaxial extensional flow at critical stages, such as entering pores or passing through a die or spinneret. Understanding the extensional rheology of viscoelastic suspensions is crucial for optimising the efficiency of such industrial processes, but remains less explored compared with shear rheology (Shaqfeh Reference Shaqfeh2019).
Previous theoretical and experimental studies have characterised the rheology of particle-free viscoelastic fluids, particularly the rapid increase in extensional viscosity due to the polymer’s coil-stretch transition (De Gennes Reference De Gennes1974). The addition of particles can influence this mechanism, leading to distinct rheological behaviours in suspensions. Recently, Jain, Einarsson & Shaqfeh (Reference Jain, Einarsson and Shaqfeh2019) examined short-term transient rheology within a narrow range of polymer relaxation times (
$\lambda$
) and concentrations (
$c$
), and Sharma & Koch (Reference Sharma and Koch2023b
) looked at steady-state rheology across a broader range of
$\lambda$
but small
$c$
. This study extends the understanding of transient rheology in dilute sphere suspensions over a wide range of
$\lambda$
and
$c$
through two computational approaches, identifying novel particle–polymer interaction mechanisms and qualitatively explaining the only existing experimental dataset on this subject (Hall et al. Reference Hall, McKinley, Erni, Soulages and Magee2009; Soulages et al. Reference Soulages, McKinley, Hall, Magee, Chamitoff and Fincke2010; McKinley et al. Reference McKinley, Haward, Jaishankar and Soulages2011; Jaishankar et al. Reference Jaishankar, Haward, Hall, Magee and McKinley2012; McKinley & Jaishankar Reference McKinley and Jaishankar2013).
In a Newtonian fluid, in the absence of inertia, the particle stresslet (Batchelor Reference Batchelor1970) leads to additional stress or viscosity in the suspension (Einstein Reference Einstein1905; Batchelor Reference Batchelor1970; Batchelor & Green Reference Batchelor and Green1972). The stresslet represents the resistance of a rigid particle to deformation by the fluid and depends on the additional traction created on the particle surface by the fluid. In viscoelastic fluids, which consist of polymers dissolved in a Newtonian solvent, the stresslet also includes the effect of polymeric stress on the particle surface. The dissolved polymers are stretched differently in the presence of particles due to velocity perturbations. This accumulated effect of extra polymer stress throughout the fluid leads to additional stress in the suspension, termed particle-induced polymer stress (PIPS). The changing polymer configuration locally alters the solvent stress by modifying the local pressure and velocity. The sum of the stresslet due to polymer stress and the solvent stress modification is termed the interaction stresslet. Therefore, the two-way interaction between particles and polymers creates two additional stress mechanisms in the suspension: PIPS, where particles modify the polymer stretch or configuration, and the interaction stresslet, where polymers exert additional stress on the particle. The net particle–polymer interaction stress is the sum of PIPS and the interaction stresslet.
The stress measured in an experiment with a homogeneous suspension is an average over all possible particle configurations (Hinch Reference Hinch1977). The stresslet, or the average of the extra stress inside the rigid particles, can be expressed in terms of a surface integral yielding the first moment of the fluid stress over the particle surface, as shown by Batchelor (Reference Batchelor1970). This approach renders the knowledge of the constitutive relation within the particle unnecessary. In sufficiently dilute suspensions, particle–particle interactions are rare. Therefore, the flow and polymer configuration field in an unbounded fluid around an isolated particle fully determine the rheology of the suspension. While PIPS represents the extra polymeric stress in the fluid, previous attempts by Greco, D’Avino & Maffettone (Reference Greco, D’Avino and Maffettone2007), Housiadas & Tanner (Reference Housiadas and Tanner2009) and Jain et al. (Reference Jain, Einarsson and Shaqfeh2019) to approximate PIPS through volume averaging have been inappropriate. As discussed by Koch, Lee & Mustafa (Reference Koch, Lee and Mustafa2016), based on the far-field Newtonian velocity and strain rate disturbance scalings of
$1/r^2$
and
$1/r^3$
, respectively, the additional polymeric stress in a steady linear flow around a sphere decays as
$1/r^3$
at large distances
$r$
from the particle. Consequently, volume averaging leads to logarithmic divergence, which causes slow or absent convergence in direct numerical simulations. Koch et al. (Reference Koch, Lee and Mustafa2016) showed that this logarithmic divergence arises from the polymer stress linearised about the velocity and polymer stress in the fluid undisturbed by the particle. They also demonstrated that the ensemble average of this linearised stress is zero. Therefore, removing the linearised stress from the extra stress before approximating the ensemble averaging with volume averaging leads to convergent integrals. This approach is also numerically beneficial, as it accelerates convergence with domain size by removing the slower-decaying component of the stress from the volume integral. In a transient flow, the initial extra polymer stress is small, and equivalent PIPS is obtained with or without removing the linearised stress. However, as the extra polymer stress increases in both magnitude and spatial extent, the divergence in the linearised stress appears at a time proportional to the polymer relaxation time.
In numerical or theoretical studies of viscoelastic fluids, polymers are typically modelled using a continuum equation for the polymer configuration tensor, which is proportional to the polymer stress (Bird & Giacomin Reference Bird and Giacomin2016). A simple way to capture polymer elasticity is through the dumbbell model, which represents a linear polymer molecule as a spring with two Brownian beads influenced by velocity gradients. The equations for the polymer configuration tensor,
$\boldsymbol{\varLambda}$
, are derived from the outer product of the dumbbell’s end-to-end vector with itself and then averaged over all possible polymer orientations to obtain continuum-level governing equations. A nonlinearly elastic spring is used to capture the finite extensibility of the polymer molecule, which is crucial in strong flows such as uniaxial extension. In simpler models like the Oldroyd-B model, the polymer undergoing extensional flow stretches indefinitely, and the polymer stress can become unbounded. The FENE (finite extensible nonlinear elastic) models, particularly the FENE-P model (Peterlin Reference Peterlin1968), effectively capture the coil-stretch transition at small polymer concentrations,
$c$
(Anna & McKinley Reference Anna and McKinley2008), which is a key mechanism in this study. For higher
$c$
, the Giesekus model is more appropriate as it better describes the behaviour of high-concentration viscoelastic liquids (Khan & Larson Reference Khan and Larson1987).
In the steady state of a liquid with dilute polymer concentration,
$c$
, particles and polymers interact differently before and after the coil-stretch transition of the undisturbed polymers (Sharma & Koch Reference Sharma and Koch2023b
), which occurs at
$De = 0.5$
. The Deborah number,
$\textit{De}$
, is defined as the product of polymer relaxation time,
$\lambda$
, and imposed extension rate. When undisturbed polymers are in a coiled state (
$De \lt 0.5$
), velocity perturbations created by the particle cause extra stretching (beyond that caused by the imposed flow) around the extensional axis downstream of the particle, leading to a wake of highly stretched polymers. This results in additional stress due to particle–polymer interaction. When polymers are already in their stretched state, these regions of extra stretching do not significantly alter the polymer stretch near the particle compared with the far field. However, low-stretching regions near the stagnation points on the particle surface, at the compressional (front) and extensional (rear) axes, cause the polymers to collapse nearly to their equilibrium or unstressed state. This creates regions of collapsed polymers around the particle surface. Therefore, beyond the coil-stretch transition of undisturbed polymers (
$De \gt 0.5$
), even at low
$c$
, the significant particle–polymer interaction stress can lead to a reduction in the suspension viscosity compared with the viscosity of the particle-free polymeric fluid.
Recent numerical results by Jain et al. (Reference Jain, Einarsson and Shaqfeh2019) have explored the extensional rheology of dilute particle-filled viscoelastic liquids. Our study enhances these results in several ways. Firstly, while they used volume averaging to obtain the suspension stress equivalent to PIPS, we provide a quantitative correction using the more appropriate formulation by Koch et al. (Reference Koch, Lee and Mustafa2016). Secondly, we examine a broader range of parameters to characterise the suspension rheology. Although Jain et al. (Reference Jain, Einarsson and Shaqfeh2019) provided results for three different polymer constitutive models, their simulations were limited in terms of maximum time or Hencky strain (time non-dimensionalised with strain rate),
$H$
, and the range of
$\textit{De}$
. Their simulations only extended up to the coil-stretch transition of undisturbed polymers, missing some profound features such as the negative particle–polymer interaction stress. We observe that significant particle–polymer interactions continue beyond the coil-stretch transition. In addition to considering larger
$H$
, we explore important qualitative and quantitative variations in suspension rheology by examining a wider range of
$\textit{De}$
,
$c$
, maximum polymer extensibility (
$L$
) in FENE-P and polymer mobility (
$1/\alpha$
) in Giesekus liquids. Lastly, we provide more detailed mechanistic explanations underlying the rheological observations.
To fully elucidate the rheology of suspensions, we implement two synergistic methods. Firstly, we adopt a semi-analytical approach, previously utilised by Koch et al. (Reference Koch, Lee and Mustafa2016) and Sharma & Koch (Reference Sharma and Koch2023b
), which employs a regular perturbation expansion in
$c$
and a generalised reciprocal theorem. This technique assesses the rheological behaviour of particle suspensions in viscoelastic liquids with low
$c$
values, across various
$\textit{De}$
, maximum polymer extensibility (
$L$
), and
$H$
. Its computational efficiency is due to the straightforward numerical integration of ordinary differential equations obtained via the method of characteristics, allowing for exceptional spatial and temporal resolution and the capability to handle large computational domains. From this approach, we derive critical scaling relationships in
$L$
and identify key mechanisms influencing suspension rheology across a wide range of
$\textit{De}$
. Building on these insights, we then employ direct numerical simulations to investigate the rheology at higher
$c$
values, where novel particle–polymer interaction mechanisms emerge. These simulations are facilitated by our in-house finite difference solver, introduced in Sharma & Koch (Reference Sharma and Koch2023a
). This solver, formulated in prolate spheroidal coordinates, precisely models the particle surface. The chosen coordinate system naturally concentrates the grid near the particle surface in Euclidean space, achieving high resolution in areas with significant gradients of flow variables. The computational domain’s outer boundary is defined by a far-field spherical surface where the extensional flow conditions are applied, simulating an unbounded fluid surrounding a particle.
The rest of the paper is organised as follows. Section 2 discusses the equations governing fluid velocity and polymer stress around the particle, as well as the necessary rheological quantities. Since the additional stresses in the suspension arising from particle–polymer interaction, discussed in §§ 4 and 5, depend on the coil-stretch transition of the polymers far from the particles, the evolution of the undisturbed polymer stress is briefly reviewed in § 3. This discussion also helps compare the suspension rheology across different polymer constitutive models (FENE-P and Giesekus) used in our study. Section 4 characterises the transient nature of the particle–polymer interaction on suspension rheology for dilute
$c$
using a semi-analytical method. The effects of polymer concentration on suspension rheology using direct numerical simulations are illustrated in § 5. We summarise the only available experimental results on the effect of particles on the extensional flow of polymeric fluids in § 6 and relate these to our computational findings. The main conclusions are provided in § 7. More details about the rheological quantities, interaction stresslet, and PIPS than mentioned in § 2 are provided in Appendix A, where the treatment of ensemble averaging by first removing the linearised stress is also discussed. Appendices B and C provide details about the semi-analytical method and the techniques and validation for the direct numerical simulations, respectively.
2. Mathematical formulation and different stresses arising in a particle suspension in viscoelastic liquids
We consider the inertia-less and incompressible polymeric fluid flow in a dilute suspension of spheres. Thus, throughout the suspension, mass and momentum conservation imply a divergence-free velocity,
$\boldsymbol{u}$
, and stress,
$\boldsymbol{\sigma}$
, field, respectively,
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0,\hspace {0.2in}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\sigma}=0. \end{equation}
Since the fluid consists of polymers dissolved in a Newtonian solvent, the stress is the sum of a Newtonian solvent stress,
$\boldsymbol{\tau}$
, and a polymeric/viscoelastic stress,
$\boldsymbol{\varPi}$
,
\begin{equation} \boldsymbol{\sigma}=\boldsymbol{\tau}+c\boldsymbol{\varPi}=-p\boldsymbol{\delta}+2\boldsymbol{e}+c\boldsymbol{\varPi}. \end{equation}
Here,
$p$
is the hydrodynamic pressure,
$\boldsymbol{e} = (\boldsymbol{\nabla}\boldsymbol{u} + (\boldsymbol{\nabla}\boldsymbol{u})^{{T}})/2$
is the strain rate tensor and
$\boldsymbol{\delta}$
is the identity tensor. According to the Huggins equation, for dilute polymer solutions, the volume or mass concentration of the polymers is linearly proportional to the ratio of polymer to solvent viscosity (Rubinstein & Colby Reference Rubinstein and Colby2003). In non-dilute solutions, while remaining a monotonic function, the viscosity ratio becomes a quadratic function of the volume concentration. In this work, to maintain simplicity, we refer to the ratio of polymer to solvent viscosity as polymer concentration,
$c$
. We consider the effects of suspension under uniaxial extensional flow, such that the boundary conditions on (2.1) are
\begin{equation} \boldsymbol{u}=0,\text{ on particle surface, and, }\boldsymbol{u}=\langle \boldsymbol{u}\rangle \text{ as} |\boldsymbol{r}|\rightarrow \infty , \end{equation}
where the imposed velocity is averaged over an ensemble of all particle configurations with this average indicated by the angular brackets
\begin{equation} \langle \boldsymbol{u}\rangle =\boldsymbol{r}\boldsymbol{\cdot}\langle \boldsymbol{e}\rangle ,\hspace {0.2in} \langle {e_{\textit{ij}}}\rangle =\delta _{i1}\delta _{j1}-\tfrac{1}{2}(\delta _{i2}\delta _{j2}+\delta _{i3}\delta _{j3}). \end{equation}
The radius of the spheres, the inverse of the imposed extensional rate,
$\dot {\epsilon}$
, and the solvent viscosity are used as the length, time and viscosity scales to non-dimensionalise the governing equations.
We model the polymeric stress with either FENE-P or Giesekus constitutive equations, which consider the polymer to be a dumbbell as mentioned in § 1. In a uniaxial extensional flow, the FENE-P model is appropriate for capturing the polymer stresses in liquids with dilute polymer concentrations (Anna & McKinley Reference Anna and McKinley2008), whereas the Giesekus model is useful for concentrated solutions (Khan & Larson Reference Khan and Larson1987). The polymer’s stress is a function of its conformation,
$\boldsymbol{\varLambda}=\langle \boldsymbol{q}\boldsymbol{q}\rangle _{\textit{polymer orientation}}$
, defined as the average over orientations of the outer product of the relative vector of the ends of the dumbbell,
$\boldsymbol{q}$
, with itself. The stress-conformation relation for each model is
\begin{equation} \boldsymbol{\varPi}=\begin{cases} \dfrac {1}{De}(f\boldsymbol{\varLambda}-b\boldsymbol{\delta}), f=\dfrac {L^2}{L^2-\text{tr}(\boldsymbol{\varLambda})}, \hspace {0.2in} b=\dfrac {L^2}{L^2-\text{tr}(\boldsymbol{\delta})},&\text{FENE-P},\\[8pt] \dfrac {1}{De}(\boldsymbol{\varLambda}-\boldsymbol{\delta}),&\text{Giesekus}. \end{cases} \end{equation}
The underlying velocity field influences
$\boldsymbol{\varLambda}$
according to the following equation:
\begin{equation} \frac {\partial \boldsymbol{\varLambda}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\varLambda}=\boldsymbol{\nabla}\boldsymbol{u}^{{T}}\boldsymbol{\cdot}\boldsymbol{\varLambda}+\boldsymbol{\varLambda}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}+\begin{cases} \dfrac {1}{De}(b\boldsymbol{\delta}-f\boldsymbol{\varLambda}),&\text{FENE-P}\\[8pt] \dfrac {1}{De}((\boldsymbol{\delta}-\boldsymbol{\varLambda})-\alpha (\boldsymbol{\delta}-\boldsymbol{\varLambda})\boldsymbol{\cdot}(\boldsymbol{\delta}-\boldsymbol{\varLambda})),&\text{Giesekus}. \end{cases} \end{equation}
In the FENE-P model,
$L$
is the maximum extensibility of the polymer normalised with its radius of gyration. In the Giesekus model,
$0\leqslant \alpha \leqslant 1$
is the mobility parameter that models the ability of the entangled polymers to restrict the motion and stress of one another. These dumbbell models capture the polymer convection by the underlying velocity field (second term on the left), its stretching and rotation by the underlying velocity gradient (first two terms on the right), and its relaxation to the equilibrium state,
$\boldsymbol{\varLambda}_{{equil.}}=\boldsymbol{\delta}$
(last term on the right, modelled differently in FENE-P and Giesekus).
In the dumbbell models, the polymer stretch,
$\mathcal{S}$
, or the 2-norm of the relative vector of the dumbbell’s ends is
\begin{equation} \mathcal{S}= \langle ||\boldsymbol{q}||_2\rangle _{\textit{polymer orientation}}=\sqrt {\text{tr}(\boldsymbol{\varLambda})}. \end{equation}
A non-dimensional group, the Deborah number (
$\textit{De}$
), appears in these equations and is the ratio of polymer relaxation time,
$\lambda$
, to the flow characteristic time,
$1/\dot {\epsilon}$
,
\begin{equation} De=\lambda \dot {\epsilon}. \end{equation}
Neglecting fluid inertia leads to quasi-steady velocity and pressure fields that respond instantaneously to the polymer stress in (2.1). This requires the momentum diffusion time through the fluid (with kinematic viscosity
$ \nu$
) in the region surrounding a particle of radius
$ a$
, given by
$ a^2/\nu$
, to be small compared with
$ \lambda$
; i.e. the elasticity number (the ratio of Deborah to Reynolds number),
$El = De/Re = \lambda \nu / a^2 \gg 1$
, must be large – a common limit in polymeric fluids.
While (2.2) captures the fluid stress, in a suspension of particles additional stress components arise whose physical origins are worth noting. The stress in a suspension of particles is the average stress over the ensemble of all possible particle configurations. In a particular instance of particle configuration, the stress at any point within the suspension is the sum of the stress given by (2.2) and the extra stress within the particle phase (which is zero in the fluid phase). In the rheology of an incompressible suspension, the deviatoric or traceless part of the suspension stress is most relevant, as the trace can be absorbed into a modified pressure. The ensemble-averaged deviatoric suspension stress in a particle suspension within a viscoelastic liquid is
\begin{equation} \langle \boldsymbol{\hat{\sigma}}\rangle =\langle {\boldsymbol{\sigma}}\rangle -\frac {1}{3}\boldsymbol{\delta}\text{tr}(\langle \boldsymbol{\sigma}\rangle )=2\langle \boldsymbol{e}\rangle +c\langle \boldsymbol{\hat{\varPi}}\rangle +n\boldsymbol{\hat{S}}(\boldsymbol{\sigma}), \end{equation}
or
\begin{equation} \langle \boldsymbol{\hat{\sigma}}\rangle =(2+5\phi )\langle \boldsymbol{e}\rangle +c\boldsymbol{\hat{\varPi}}{}^{U}+{c\phi}(\boldsymbol{\hat{\varPi}}{}^{\textit{PP}}+\boldsymbol{\hat{S}}{}^{\textit{PP}}). \end{equation}
Since
$\boldsymbol{\hat{\sigma}}$
represents the stress field at any point within the suspension, its ensemble average incorporates the effects from both the particle and fluid phases. Here,
$n$
is the number of particles per unit volume and
$\phi =nV_p/(V_f+nV_p)$
is the particle volume fraction where the volume of each particle is
$V_p$
and the total suspension volume is the sum of particle (
$nV_p$
) and fluid volume (
$V_f$
). The term
$2\langle \boldsymbol{e}\rangle$
represents the imposed rate of strain,
$c\langle \boldsymbol{\hat{\varPi}}\rangle =c(\boldsymbol{\hat{\varPi}}{}^{U}+\phi \boldsymbol{\hat{\varPi}}{}^{\textit{PP}})$
is the ensemble average of the polymeric stress in the fluid and
$n\boldsymbol{\hat{S}}(\boldsymbol{\sigma})=5\phi \langle \boldsymbol{e}\rangle +c\phi \boldsymbol{\hat{S}}{}^{\textit{PP}}$
arises due to the extra stress in the particle phase and is termed as the particle stresslet. The third term in (2.10) is the polymer stress in the absence of the particles, and the second, fourth and fifth terms in (2.10) represent the effect of particles. The second term,
$5\phi \langle \boldsymbol{e}\rangle$
, is the Newtonian stresslet in a suspension of spheres (first calculated by Einstein (Reference Einstein1905)). The last two terms,
${c\phi}(\boldsymbol{\hat{\varPi}}{}^{\textit{PP}}+\boldsymbol{\hat{S}}{}^{\textit{PP}})$
, represent the particle–polymer interaction stress in the suspension. The presence of particles leads to a deviation in the fluid velocity that alters the polymer stress within the fluid, and the extra stress hence created manifests as the PIPS,
${c\phi}\boldsymbol{\hat{\varPi}}{}^{\textit{PP}}$
. The presence of polymers creates an additional stress in the particle phase which is represented by the interaction stresslet,
${c\phi}\boldsymbol{\hat{S}}{}^{\textit{PP}}$
.
Due to the symmetry of the imposed uniaxial extensional flow around a sphere, we may express
\begin{equation} \langle \boldsymbol{\hat{\sigma}}\rangle =2\mu _{\textit{ext}}\langle \boldsymbol{e}\rangle , \text{with, }\mu _{\textit{ext}}= 1+2.5\phi +0.5 c \hat{\varPi}{}^{U}\bigg[1+{\phi}\frac {\hat{\varPi}^{\textit{PP}}+\boldsymbol{\hat{S}}{}^{\textit{PP}}}{\hat{\varPi}{}^{U}}\bigg], \end{equation}
being the extensional viscosity of the suspension. It is the sum of the extensional viscosity of the polymeric fluid,
$\mu _{\textit{fluid}}$
, and the net effect of the addition of particles to the suspension stress,
$\mu _{\textit{part}}$
\begin{equation} \mu _{\textit{fluid}}=1+0.5c\hat{\varPi}{}^{U},\hspace {0.2in}\mu _{\textit{part}}=\phi (2.5+0.5c(\hat{\varPi}^{\textit{PP}}+\boldsymbol{\hat{S}}{}^{\textit{PP}})). \end{equation}
Further details on the constituents of particle–polymer interaction stress,
${c\phi}(\boldsymbol{\hat{\varPi}}{}^{\textit{PP}}+\boldsymbol{\hat{S}}{}^{\textit{PP}})={c\phi}(\hat{\varPi}^{\textit{PP}}+\boldsymbol{\hat{S}}{}^{\textit{PP}})\langle \boldsymbol{e}\rangle$
are provided in Appendix A. We will consider the extensional component of the stress when illustrating the polymer stress for both the undisturbed fluid,
$\boldsymbol{\hat{\varPi}}{}^{U}$
, in figures 1 and 2 in § 3 and that due to the effect of particle polymer interaction
$\hat{\varPi}^{\textit{PP}}+\boldsymbol{\hat{S}}{}^{\textit{PP}}$
in figures 3–5, and 11–18 in §§ 4 and 5.
3. Undisturbed polymer stress and the coil-stretch transition phenomenon
In extensional rheology of a startup flow with constant strain rate, time is non-dimensionalised with the imposed extension rate and is termed as Hencky strain,
$H$
. In a uniaxial extensional flow, the transient behaviour of the undisturbed polymer stress,
$c\hat{\varPi}{}^{U}\langle \boldsymbol{e}\rangle$
, is significant and qualitatively impacts the particle–polymer interaction stress. As we will later see, the polymer behaviour in various regions around the particle relative to the undisturbed polymers is crucial for understanding the particle–polymer interaction stress. The undisturbed polymer stress varies linearly with
$c$
for the entire range of relevant parameters, so we consider this stress normalised with
$c$
, i.e.
$\hat{\varPi}{}^{U}\langle \boldsymbol{e}\rangle$
, for the FENE-P and Giesekus liquids in this section.
Since the undisturbed stress is spatially homogeneous for imposed linear flows, such as the uniaxial extensional flow considered here, the governing equations for this stress are ordinary differential equations (ODEs). Specifically, the convective term in (2.6) is absent for the undisturbed polymer configuration equations. We solve the corresponding ODEs using MATLAB’s ode15s for the small
$c$
results in § 4 and a variable time step second-order Runge–Kutta method for the results at an arbitrary
$c$
in § 5.
3.1. The FENE-P model
The evolution of undisturbed polymer stress in a FENE-P polymeric liquid is shown in figure 1 for different Deborah numbers,
$\textit{De}$
, and maximum polymer extensibility,
$L$
. Initially, the polymers are in equilibrium (zero stress) and are subsequently stretched due to the imposed extension, causing
$\hat{\varPi}{}^{U}$
to increase monotonically with
$H$
. For
$De\lt 0.5$
, the polymer stretch,
$\sqrt {\text{tr}(\boldsymbol{\varLambda}^{U})}$
, is much smaller than
$L$
, resulting in almost identical curves for different
$L\geqslant 50$
at
$De=0.2$
and 0.4 in figure 1. Additionally, because
$\sqrt {\text{tr}(\boldsymbol{\varLambda}^{U})}\ll L$
in the
$De\lt 0.5$
regime for all
$H$
, increasing
$\textit{De}$
causes the steady state to be reached later, as larger
$\textit{De}$
implies a larger final polymer stretch.
Figure 1. Evolution of undisturbed (no particles) polymer stress,
$\hat{\varPi}{}^{U}$
with Hencky strain,
$H$
for: (a)
$De=0.2$
, (b)
$De=0.4$
, (c)
$De=0.6$
, (d)
$De=2.0$
and (e)
$De=5.0$
for different
$L$
. For the latter three cases
$\hat{\varPi}{}^{U}$
is normalised with
$L^2$
and
$H$
with
$\log (L)$
. All cases share the same legend for
$L$
shown in (a).
For
$De\geqslant 0.5$
, as shown in figure 1, we observe the well-known coil-stretch transition (De Gennes Reference De Gennes1974; Perkins, Smith & Chu Reference Perkins, Smith and Chu1997) – the polymer stress increases much more rapidly with
$H$
beyond a certain point as the coiled polymer chains uncoil and approach full extension. Consequently, the polymer stress (proportional to the polymer configuration, which is the outer product of the polymer end-to-end vector) increases as
$L^2$
. Therefore, the
$L^2$
-normalised graphs for sufficiently large
$L$
at
$De=0.6$
, 2.0 and 5.0 are identical in the steady state. The polymers stretch towards
$L$
at a rate inversely proportional to the logarithm of
$L$
. The transient behaviour of different curves is similar for a normalised Hencky strain,
$H/\log (L)$
, for very large
$\textit{De}$
, as shown here for 2.0 and 5.0 at
$L\geqslant 50$
. In this regime of large
$\textit{De}$
, polymers stretch faster towards the steady state as
$\textit{De}$
increases, since the polymers have a maximum limit,
$L$
, that can be approached more quickly with a larger extension rate.
The effects of limited polymer extensibility are observed close to the coil-stretch transition point of
$De=0.5$
. For
$De=0.4$
, with
$\sqrt {\text{tr}(\boldsymbol{\varLambda}^{U})}\ll L$
, the collapsed polymers stretch much less for
$L=10$
than larger
$L$
, resulting in a lower steady-state
$\hat{\varPi}{}^{U}$
. For
$De=0.6$
, polymers undergo the coil-stretch transition, and
$\text{tr}(\boldsymbol{\varLambda}^{U})\sim L^2$
, but the stretch of larger
$L$
polymers is still not as close to their respective
$L$
as it is for
$L=10$
. When
$\textit{De}$
is larger (
$De=$
2.0 and 5.0 in figure 1), the coil-stretch transition is more rapid. Due to the large extension rate, even the larger
$L$
polymers get extended very close to their respective
$L$
, and the (normalised) steady state at all
$L$
is similar.
3.2. The Giesekus model
The evolution of the undisturbed polymer stress for different levels of polymer entanglement (proportional to the mobility parameter
$\alpha$
) predicted by the Giesekus model is shown as solid lines for three different
$De=0.4$
, 2.0, and 5.0 in figure 2. While the FENE-P model imposes a maximum polymer extensibility explicitly, the Giesekus model implicitly restricts the polymer’s maximum extension through the mobility parameter
$\alpha$
, which increases for more entangled polymers. By analysing the relaxation term in (2.6), we can define an effective maximum extensibility,
$L_{\textit{Giesekus}}=\alpha ^{-0.5}$
. This equivalence between the FENE-P and Giesekus models was previously known. Housiadas & Beris (Reference Housiadas and Beris2013) obtained it through the expression for the steady-state extensional viscosity at large
$\textit{De}$
. Therefore, in figure 2, we also show the curves obtained from the FENE-P model as dashed lines with
$L=\alpha ^{-0.5}$
for the same values of
$\alpha$
as those considered in the solid lines representing the Giesekus model. The qualitative features of the undisturbed stress in Giesekus fluids are similar to those in FENE-P fluids discussed in § 3.1 above. The stress is independent of
$\alpha$
(or
$L_{\textit{Giesekus}}$
) for
$De\lt 0.5$
for small enough
$\alpha$
(large enough
$L_{\textit{Giesekus}}$
) throughout the transient and scales as
$\alpha ^{-1}$
or
$L_{\textit{Giesekus}}^2$
above
$De=0.5$
in the steady state. In the large
$\textit{De}$
regime, the polymers stretch towards
$L$
at a rate inversely proportional to the logarithm of
$\alpha$
. The steady-state polymer stress of a Giesekus fluid for small enough
$\alpha$
(large enough
$L$
) in the small
$\textit{De}$
regime and in the large
$\textit{De}$
regime for all
$\alpha$
considered is the same as that for a FENE-P fluid with
$L=\alpha ^{-0.5}$
. However, in the large
$\textit{De}$
regime, it takes longer for the polymers in a Giesekus fluid to reach the steady state relative to those in a FENE-P fluid. This is likely because, in addition to uncoiling, the Giesekus polymers also have to disentangle from one another.
Figure 2. Evolution of undisturbed (no particles) polymer stress,
$\hat{\varPi}{}^{U}$
with Hencky strain,
$H$
for: (a)
$De=0.4$
, (b)
$De=2.0$
and (c)
$De=5.0$
for different
$\alpha$
in the Giesekus model (solid lines) and equivalent
$L=\alpha ^{-0.5}$
in the FENE-P model (dashed lines). For (b) and (c)
$\hat{\varPi}{}^{U}$
is normalised with
$\alpha ^{-1}$
and
$H$
with
$-0.5 \log (\alpha )$
. All cases share the same legend for
$\alpha$
or
$L$
shown in (a).
4. Suspension rheology of dilute polymeric solutions (FENE-P liquids)
At very low polymer concentrations (
$c$
), the polymer configuration is primarily influenced by the velocity fluctuations generated by the particles in a Newtonian fluid. This means that the Newtonian velocity field is sufficient to describe the polymer behaviour. This simplification allows for more efficient calculations than those using computational fluid dynamics (CFD) techniques, which involve discretisation of the equations governing fluid mass and momentum along with the polymer constitutive equations. Here, the extensional viscosity of the suspension, denoted as
$\mu _{\textit{ext}}$
in (2.11), is obtained using a combination of regular perturbation in
$c$
, a generalised reciprocal theorem, and a method of characteristics.
The leading-order polymer constitutive equation is obtained by solving the evolution of polymers along the streamlines (characteristics) of the Newtonian flow. The leading-order effect of polymers on the fluid velocity and pressure field, which modifies the particle stresslet, is evaluated using a generalised reciprocal theorem. This approach circumvents the need for a numerical solution of the
$\mathcal{O}(c)$
mass and momentum equations. The use of the method of characteristics makes this method similar to Lagrangian particle methods (LPM) for tracer transport (Bosler et al. Reference Bosler, Kent, Krasny and Jablonowski2017), in contrast to CFD methods used more commonly in Eulerian reference frame. However, unlike the use of LPM in a general flow scenario, here the fluid velocity is steady, making the implementation more straightforward. This semi-analytical methodology is similar to that considered by Koch et al. (Reference Koch, Lee and Mustafa2016) and Sharma & Koch (Reference Sharma and Koch2023b
). Further details and the required changes due to the transient nature of polymer evolution in a uniaxial extensional flow are discussed in Appendix B.
The results of these calculations, spanning a wide range of Deborah numbers (
$\textit{De}$
) and maximum polymer extensibility (
$L$
), are presented in § 4.1 before the underlying mechanisms are elucidated in § 4.2.
4.1. Rheological observations
The extensional viscosity of the suspension and its contributions from the particle–polymer stress discussed in this section are formally
$\mathcal{O}(c)$
and are calculated from (B2) to (B4). Since we only consider
$\mathcal{O}(c)$
effects, for brevity of notation, we do not include
$0$
in the superscript of various quantities appearing in these equations while discussing the following results.
Figure 3. Evolution of total particle–polymer interaction stress,
$(\hat{\varPi}^{\textit{PP}}+\boldsymbol{\hat{S}}{}^{\textit{PP}})$
with Hencky strain,
$H$
for: (a)
$De=0.4$
, (b)
$De=0.6$
and (c)
$De=5.0$
. For the latter two cases
$H$
is normalised with
$\log (L)$
and stress with
$L^2$
. All cases share the same legend for
$L$
shown in (a).
4.1.1. Net effect of particle–polymer interaction on the extensional viscosity of the suspension
figure 3 illustrates the evolution of the net particle–polymer interaction stress,
$\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}$
, for different Deborah numbers (
$De = 0.4$
, 0.6 and 5.0) and maximum polymer extensibilities (
$L = 10$
, 50, 100 and 200) with respect to the Hencky strain,
$H$
. In the case of
$De = 0.6$
and 5.0,
$H$
is normalised with
$\log (L)$
, and the stress is normalised with
$L^2$
. To focus specifically on the interaction stress, figure 4 presents the evolution of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
, i.e. the particle–polymer interaction stress normalised with
$\hat{\varPi}{}^{U}$
from figure 1. This allows for a clearer understanding of the relative contribution of the interaction stress compared with the undisturbed stress. The figure includes data for
$De = 0.2, 0.4, 0.6, 2.0$
and 5.0, and for each Deborah number, it shows the results for
$L = 10$
, 50, 100 and 200. At
$H = 0$
,
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U} = 2.5$
, consistent with the small time estimates obtained by (A16). For
$De = 0.2$
and 0.4, the stress at large
$L \geqslant 50$
is independent of
$L$
for all values of
$H$
, while for
$De = 0.6$
, 2.0 and 5.0, the stress scales with
$L^2$
for all values of
$H / \log (L)$
. The attainment of steady state occurs more slowly for larger Deborah numbers below
$De = 0.5$
, but faster for
$De \gt 0.5$
. These scaling relationships and trends with respect to
$L$
and
$H$
observed in the particle–polymer interaction stress are similar to those observed in the undisturbed polymer stress shown in figure 1 and discussed in § 3.
For
$De=0.2$
we observe a monotonic increase in the normalised interaction stress at each
$L$
in figure 4(a). This indicates that particle–polymer interaction contributes to an increase in the extensional viscosity,
$\mu _{\textit{ext}}$
, of the suspension with the Hencky strain until reaching a steady state. This positive effect arises from the net contribution of particles,
$\mu _{\textit{part}}$
(2.12). At the lowest
$L=10$
shown, the normalised interaction stress behaves similarly to the higher
$L$
cases but with a lower magnitude. Similar to the undisturbed stress,
$\hat{\varPi}{}^{U}$
, in figure 1, the normalised interaction stress,
$(\hat{\varPi}^{\textit{PP}}+\boldsymbol{\hat{S}}{}^{\textit{PP}})/\hat{\varPi}{}^{U}$
in figure 4 is larger for
$De=0.4$
than
$De=0.2$
. From this observation and additional data for
$De = 0.1$
and 0.3 (not shown), we conclude that in the
$De \lt 0.5$
regime, increasing
$\textit{De}$
enhances the positive contribution of
$\mu _{\textit{part}}$
to
$\mu _{\textit{ext}}$
at all
$H$
. For
$De = 0.4$
, the interaction of particles with
$L = 10$
polymers is qualitatively different from that with higher
$L$
polymers, as shown by the transient peak for
$L = 10$
in figure 4(b).
Figure 4. Evolution of total interaction stress normalised with the undisturbed (no particles) polymer stress,
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
, with Hencky strain,
$H$
, for: (a)
$De = 0.2$
, (b)
$De = 0.4$
, (c)
$De = 0.6$
, (d)
$De = 2.0$
and (e)
$De = 5.0$
for different
$L$
. For the latter three cases,
$H$
is normalised with
$\log (L)$
. All cases share the same legend for
$L$
shown in (a).
The effect of particle–polymer interaction becomes more profound and interesting for
$\textit{De}$
beyond the undisturbed coil-stretch transition value of 0.5. In this regime, the evolution of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
with
$H$
is non-monotonic, as shown in the
$De = 0.6$
, 2.0 and 5.0 plots in figure 4. Initially, at smaller
$H$
, it increases rapidly to large positive values and then sharply falls to negative values. Before settling into a negative steady state for all
$L$
at
$De = 2.0$
and 5.0, there is a valley in
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
(and also in
$\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}$
as depicted in figure 3(c)). The magnitude of this valley increases with
$\textit{De}$
. The sharp drop to negative values occurs around the same
$H$
as when undisturbed polymers undergo the coil-stretch transition (figure 1).
We recently considered the steady-state extensional rheology in Sharma & Koch (Reference Sharma and Koch2023b
). In that study, we found that the interaction stress
$\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}$
increases with
$\textit{De}$
in the
$De \lt 0.5$
regime and is independent of
$L$
for
$L \geqslant 50$
. However, for a value of
$\textit{De}$
in the range
$0.55 \lessapprox De \lessapprox 0.65$
, depending on
$L$
,
$\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}$
changes from positive to negative and scales as
$L^2$
for
$L \geqslant 50$
(for larger
$\textit{De}$
, the
$L^2$
scaling extends to a lower
$L = 10$
). Following this negative value, further increases in
$\textit{De}$
make
$\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}$
more negative. In the current study, the same conclusions about the steady-state effect of the particle–polymer interaction as Sharma & Koch (Reference Sharma and Koch2023b
) are obtained from the largest
$H$
values shown here.
From figure 3(b), we find the maximum and steady-state values of interaction stress
$\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}$
to be
$0.1L^2$
and
$-0.33L^2$
for
$De = 0.6$
and
$L = 100$
. The maximum, minimum, and steady-state values for
$De = 5.0$
are
$7.1L^2$
,
$-2.7L^2$
and
$-1.1L^2$
, as shown in figure 3(c). Therefore, in the initial part of the transient, the maximum net particle extensional viscosity is
$\mu _{\textit{part}} = (2.5 + 1000c)\phi$
and
$(2.5 + 71\,000c)\phi$
for
$De = 0.6$
and 5.0, respectively. Considering a very small value of
$c = 0.001$
, expected to be well within the range of validity of our low
$c$
methodology, this implies a maximum particle–polymer interaction led increase in
$\mu _{\textit{part}}$
of 0.4 and 28.4 times compared with the Newtonian
$\mu _{\textit{part}} = 2.5\phi$
. From similar calculations at the largest
$H$
considered here,
$\mu _{\textit{part}} = (2.5 - 3300c)\phi$
and
$(2.5 - 11\,200c)\phi$
for
$De = 0.6$
and 5.0, respectively. Even with a low
$c$
of 0.001, the net particle effect of
$\mu _{\textit{part}}$
is negative. This is a surprising and interesting result as it indicates that adding particles reduces the extensional viscosity of the suspension relative to the fluid viscosity, i.e.
$\mu _{\textit{ext}} \lt \mu _{\textit{fluid}}$
. For large
$\textit{De}$
values such as 2.0 and 5.0 shown here, due to the negative peak, the reduction in
$\mu _{\textit{ext}}$
is even larger in the transient phase just after its large initial increase. While the low
$c$
methodology implicitly assumes the polymer stress to remain smaller than the Newtonian stress, in the above-mentioned estimates, this assumption is violated. However, direct numerical simulations, which make no assumption on the smallness of
$c$
, in the forthcoming § 5 will demonstrate that the broad observations from this low
$c$
methodology are infact valid for
$c$
even larger than 0.001. Thus, the effect of particle–polymer interaction on the extensional rheology of a dilute suspension of spheres in a polymeric liquid is expected to be profound.
4.1.2. Decomposition of particle–polymer interaction stress
To obtain further insight into the behaviour of the interaction stress, its primary constituents, PIPS (
$\hat{\varPi}^{\textit{PP}}$
) and the interaction stresslet (
$\hat{S}^{\textit{PP}}$
), normalised with
$\hat{\varPi}{}^{U}$
, are shown in figure 5 for
$De = 0.4$
, 0.6 and 5.0 with
$L = 10$
, 50, 100 and 200 for each
$\textit{De}$
. figure 5(a–c) displays PIPS, while figure 5(d–f) shows the interaction stresslet. We do not show
$De = 0.2$
and 2.0 considered earlier, but the figures for these
$\textit{De}$
values are similar to
$De = 0.4$
and 5.0, respectively.
Figure 5. Evolution of the
$\hat{\Pi}^U$
-normalized particle-induced polymer stress (PIPS,
$\hat{\Pi}^\textit{PP}/\hat{\Pi}^U$
) and the interaction stresslet (
$\hat{S}^\textit{PP}/\hat{\Pi}^U$
) with Hencky strain
$H$
. The top row (a–c) shows PIPS, while the bottom row (d–f) shows the interaction stresslet. Each column corresponds to a different Deborah number: (a), (d)
$De=0.4$
; (b), (e)
$De=0.6$
; (c), (f)
$De=5.0$
. Different extensibility values
$L$
are indicated in the legend of panel (a), which applies to all subplots.
The qualitative features of the total interaction stress are almost fully captured by PIPS. This is observed by the similarity in figures 4(b) and 5(a), 4(c) and 5(b) and 4(e) and 5(c), which compare the total interaction stress with PIPS for
$De = 0.4$
, 0.6 and 5.0, respectively. However, the contribution of
$\hat{S}^{\textit{PP}}$
is quantitatively important during the transient. It is always positive and decreases with
$H$
for all
$L$
and
$\textit{De}$
, starting from 2.5 at
$H = 0$
. Beyond the undisturbed coil-stretch transition point of
$De = 0.5$
, the decrease in
$\hat{S}^{\textit{PP}} / \hat{\varPi}{}^{U}$
with
$H$
is not monotonic, as observed from the
$De = 0.6$
and 5.0 plots in figures 5(e) and 5( f). For large
$L$
(
$\geqslant 50$
) and
$De = 2.0$
and 5.0,
$\hat{S}^{\textit{PP}} / \hat{\varPi}{}^{U}$
decreases to a value below the steady state during the transient. In these cases, the steady-state
$\hat{S}^{\textit{PP}} / \hat{\varPi}{}^{U} = 1$
(see figure 5
f for example). With limited polymer extensibility,
$L = 10$
for
$De = 2.0$
(not shown) and
$L = 10$
and 50 for
$De = 5.0$
, there is a brief period during the transient where
$\hat{S}^{\textit{PP}} / \hat{\varPi}{}^{U}$
increases before decreasing again. Both types of intermediate non-monotonic transient behaviour of
$\hat{S}^{\textit{PP}} / \hat{\varPi}{}^{U}$
occur around the values of
$H$
when the undisturbed polymers (figure 1) are undergoing coil-stretch transition. This is also the period when
$\hat{\varPi}^{\textit{PP}} / \hat{\varPi}{}^{U}$
undergoes a rapid decrease to negative values after its rapid initial increase, as previously noted for
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
. Hence, at small
$\textit{De}$
, both PIPS and the interaction stresslet contribute in the same way towards the positive
$\mu _{\textit{part}}$
for all
$L$
and
$H$
. At large
$\textit{De}$
, when
$\mu _{\textit{part}} \lt 0$
at larger
$H$
, the effect of just PIPS is even more negative, and the positive interaction stresslet reduces PIPS’ influence.
According to the alternative decomposition of the interaction stresslet,
$\boldsymbol{\hat{S}}{}^{\textit{PP}}$
, into the undisturbed,
$\boldsymbol{\hat{S}}{}^{U} = \hat{\varPi}{}^{U}$
, and the volumetric,
$\boldsymbol{\hat{S}}{}^{\textit{Vol}}$
stresslet introduced in Appendix A.1, the value of
$\boldsymbol{\hat{S}}{}^{\textit{PP}} / \hat{\varPi}{}^{U}$
in figure 5(d–f) is simply
$1 + \boldsymbol{\hat{S}}{}^{\textit{Vol}} / \hat{\varPi}{}^{U}$
. Therefore, the variation in normalised interaction stresslet with
$H$
is fully explained by the variation of
$\boldsymbol{\hat{S}}{}^{\textit{Vol}} / \hat{\varPi}{}^{U}$
, and the alternative stresslet decomposition into the undisturbed and volumetric stresslet is a clear and direct way to understand the mechanism at play here. The volumetric stresslet (defined in (A8)) is an integral of a quantity proportional to the difference in the polymer stress relative to the undisturbed value. As the Hencky strain accumulates, the distribution of the polymer stress disturbance, and not the total polymer stress, throughout the fluid volume plays a role in the variation of the interaction stresslet normalised with undisturbed stress (A8). From figure 5( f), we notice that at large
$\textit{De}$
and large
$H$
,
$\boldsymbol{\hat{S}}{}^{\textit{PP}} / \hat{\varPi}{}^{U} \rightarrow 1$
, implying
$\boldsymbol{\hat{S}}{}^{\textit{Vol}} / \hat{\varPi}{}^{U} \rightarrow 0$
. Therefore, for polymers with larger relaxation time or at large extension rates, the polymer stress disturbance in the fluid volume does not affect the interaction stresslet. In the next section, we will observe that this is because, in a large volume around the sphere, polymers stretch similarly to the undisturbed ones.
4.2. Deciphering the interaction mechanism through polymer stretch around the particle
In this section, we gain insight into the changes in the extensional viscosity of the suspension,
$\mu _{\textit{ext}}$
, through the net effect of the particle,
$\mu _{\textit{part}}$
, discussed above, by observing the influence of the particle on the polymer configuration. As mentioned at the beginning of § 4, only the Newtonian velocity field affects the leading-order polymer configuration. Therefore, to understand why the particle changes the polymer configuration in any particular way, we must delve into the kinematics of the Newtonian velocity field around a sphere in uniaxial extensional flow. We take a brief detour to reintroduce a useful scalar diagnostic field for this purpose in § 4.2.1 before illustrating the particle–polymer interaction mechanism in § 4.2.2.
4.2.1. Kinematic diagnostic measures
In Sharma & Koch (Reference Sharma and Koch2023b
), we found the scalar diagnostic field termed as the fractional change in the local Deborah number field, or
$\Delta \textit{De}_{\textit{local}}$
,
\begin{equation} \Delta \textit{De}_{\textit{local}}=\sqrt {\frac {\boldsymbol{e}:\boldsymbol{e}}{\langle \boldsymbol{e}\rangle :\langle \boldsymbol{e}\rangle}}-1, \end{equation}
to be helpful in explaining the steady-state changes in polymer configuration by the particle. Here,
$\boldsymbol{e} = 0.5(\boldsymbol{\nabla}\boldsymbol{u} + (\boldsymbol{\nabla}\boldsymbol{u})^T)$
is the local rate of strain field, and
$\langle \boldsymbol{e} \rangle$
is the undisturbed or imposed rate of strain.
$\Delta \textit{De}_{\textit{local}}$
captures the local changes in stretching by the velocity gradient field compared with the far-field stretching. A positive
$\Delta \textit{De}_{\textit{local}}$
implies more stretching locally than in the far field, while
$\Delta \textit{De}_{\textit{local}} \lt 0$
implies less local stretching. In a uniaxial extensional flow,
$\langle \boldsymbol{e} \rangle : \langle \boldsymbol{e} \rangle = 1.5$
. The
$\Delta \textit{De}_{\textit{local}}$
field is shown in figure 6(a). From figure 6(a), we observe a region of largest (an intense red) local stretching compared with the far field around the extensional axis of the flow. Additionally, there is a stretching region (red region) near the particle surface at about
$45^\circ$
from the extensional axis. A red region centred around the compressional axis represents extra local stretching, as the streamlines on either side of the axis separate. There are two notable regions of intense reduction in local stretching relative to the far field (blue regions) around the front stagnation circle and rear stagnation point. The no-slip boundary condition on the particle surface implies that the velocity gradient goes to zero at these points. There is a less intense region of reduction in local stretching away from the particle surface at
$45^\circ$
from the axes.
Figure 6. (a) Fractional change in the local Deborah number field,
$\Delta \textit{De}_{\textit{local}}$
, due to a sphere in an imposed extensional flow used to locally diagnose the kinematics of the velocity field. (b) Multipole disturbances created by the sphere in a Newtonian fluid.
The Newtonian velocity field around a sphere comprises dipole (
$1/r^2$
) and octupole (
$1/r^4$
) disturbances relative to the undisturbed velocity field. These disturbances contribute to the features of
$\Delta \textit{De}_{\textit{local}}$
, with the velocity vectors for both shown in figure 6(b). It can be observed that the black arrows (representing the dipole velocity) increase in length along their direction, towards the particle, around the extensional axis. Consequently, the front of a fluid element convecting with only this velocity field will experience a larger dipole velocity magnitude than the back, resulting in stretching. In this region, the octupole velocity (orange arrows) decreases in magnitude along their direction, reducing the stretch of a fluid element. The net effect is less stretching from the velocity gradients very close to the particle surface due to the octupole. Further away where the octupole disturbances have rapidly decayed, the net effect is due to the dipole.
Around the compressional axis, the streamlines of the overall velocity field separate, causing a stretch in the direction of the extensional axis (
$z$
-axis) (this appears as the red region around the compressional axis in figure 6(a)). However, the overall velocity vectors reduce in magnitude as they approach the sphere, causing lower stretching in the direction of the compressional axis (
$r$
-axis) (this appears as the blue region around the compressional axis near the particle surface in figure 6(a)). From figure 6(b), we can observe that both these effects – increased stretching away from and reduced stretching near the particle – are enhanced due to the dipole velocity disturbance, while the octupole velocity disturbance causes an opposite effect. Around the compressional axis, both close to and away from the particle, the dipole effect dominates. The far-field stretching region does not play a dominant role in the upcoming discussion to elucidate the interaction mechanisms in the low
$c$
regime of the current § 4. However, this stretching along the compressional axis far upstream of the particle, created by the dipole disturbance of the Newtonian velocity field, will explain in § 5 the activation of another particle–polymer interaction mechanism at larger
$c$
.
The
$\Delta \textit{De}_{\textit{local}}$
field shown in figure 6(a) will be used in conjunction with a field describing the fractional change in the polymer stretch created by the particle, defined as
$\Delta \mathcal{S}$
,
\begin{equation} \Delta \mathcal{S} = \sqrt {\text{tr}(\boldsymbol{\varLambda})} / \sqrt {\text{tr}(\boldsymbol{\varLambda}^U)} - 1, \end{equation}
to understand how the polymer stretch field changes in the presence of a particle. Here,
$\sqrt {\text{tr}(\boldsymbol{\varLambda})}$
and
$\sqrt {\text{tr}(\boldsymbol{\varLambda}^U)}$
are the local and undisturbed polymer stretches, respectively. A positive
$\Delta \mathcal{S}$
indicates an increase in the local polymer stretch, while a negative
$\Delta \mathcal{S}$
implies a reduced local polymer stretch or polymer collapse. We performed a similar analysis in Sharma & Koch (Reference Sharma and Koch2023b
) but used the change in the polymer stretch field,
$\sqrt {\text{tr}(\boldsymbol{\varLambda})} - \sqrt {\text{tr}(\boldsymbol{\varLambda}^U)}$
. In the rest of this section, whenever comparative statements are made, they are relative to the current undisturbed polymers or the stretching of the imposed velocity field in the far field. Also, we remind the reader that the Hencky strain and time are the same according to our non-dimensionalisation and hence can be used interchangeably.
4.2.2. Interaction mechanism
In this section, we use
$De=0.2$
and
$0.4$
to illustrate the particle–polymer interaction mechanism in the
$De\lt 0.5$
regime, and
$De=0.6$
and
$5.0$
for the large
$\textit{De}$
regime. In Sharma & Koch (Reference Sharma and Koch2023b
), we found that the dominant contribution of
$\Delta \mathcal{S}$
along the extensional axis in the far field at steady state arises from changes in the extensional component of the configuration tensor,
$\varLambda^{\prime}_{zz}$
. Specifically,
$\varLambda^{\prime}_{zz}\sim z^{-\alpha}$
, with
$\alpha =1/De-2$
for
$0.2\lt De\lt 0.5$
,
$\alpha =4De-2$
for
$0.5\lt De\lt 1.25$
and
$\alpha =-3$
otherwise. Therefore, the length scale at which the particle influences the steady-state polymer stretch is largest in the intermediate range of
$0.2\lt De\lt 1.25$
. We will observe the signature of this steady-state region of influence of the particles on polymer stretch in the discussion of transient results below.
In the small
$\textit{De}$
regime, a wake of extra polymer stretch downstream of the particle around the extensional axis, due to
$\Delta \textit{De}_{\textit{local}}\gt 0$
, leads to a positive particle–polymer interaction stress at all times. In contrast, in the large
$\textit{De}$
regime, a region of polymer collapse around the particle is observed beyond a
$\textit{De}$
and
$L$
dependent time. This region originates from the
$\Delta \textit{De}_{\textit{local}}\lt 0$
region around the stagnation point on the particle at the compressional axis and leads to a negative particle–polymer interaction stress. The memory effects alter the interaction behaviour with
$\textit{De}$
within each regime, as described below.
A small
$De = 0.2$
is a case where the polymers react mostly to the local flow. figure 7 shows
$\Delta \mathcal{S}$
for
$De = 0.2$
and
$L = 100$
at two different
$H = 0.5$
and 2.0 (close to steady state). At
$H = 0$
, the polymers are in equilibrium everywhere; as
$H$
increases, the polymers further from the particle start to get stretched by the imposed extensional flow. At a small but finite
$H = 0.5$
for
$De = 0.2$
, the memory effects of polymer relaxation are limited such that the polymers respond only to the local velocity gradients. At this instant, there is an almost perfect positive correlation between the
$\Delta \mathcal{S}$
field in figure 7(a) and the
$\Delta \textit{De}_{\textit{local}}$
field in figure 6(a). This indicates a good local kinematic diagnostic capability of the latter. Upon further increasing
$H$
to 2.0, the polymers are stretched further in the undisturbed regions and even more so in the regions of larger local stretch. This shows itself as more intense red regions in figure 7(b) compared with figure 7(a). A more intense blue region indicates that locally in the blue or
$\Delta \textit{De}_{\textit{local}} \lt 0$
regions, the polymers are even more collapsed compared with the currently undisturbed polymers at
$H = 2.0$
than they were at
$H = 0.5$
. The weak memory effects due to a small but non-zero
$De = 0.2$
can be observed in figure 7(b), where the red or highly stretched polymer region around the extensional axis is elongated compared with the red or
$\Delta \textit{De}_{\textit{local}} \gt 0$
region in figure 6(a) around the extensional axis. This occurs because once the polymers get stretched during
$H \lt 2.0$
, they need some time or downstream distance along their convecting streamlines to relax from this stretched state. The
$\Delta \mathcal{S}$
shown in figure 7(b) has already reached the steady state for
$De = 0.2$
and
$L = 100$
.
Figure 7. Value of
$\Delta \mathcal{S}$
for
$De = 0.2$
and
$L = 100$
for Hencky strain,
$H$
: (a) 0.5 and (b) 2.0. Compared with the currently undisturbed polymers, in the red regions the polymers are stretched more, and in the blue regions they are collapsed.
The polymer stretch behaviour is qualitatively similar for
$De = 0.4$
and
$L = 100$
in figure 8 as that discussed above for
$De = 0.2$
and
$L = 100$
, but with a few differences. At
$H$
smaller than shown in figure 8,
$\Delta \mathcal{S}$
for
$De = 0.4$
is similar to the
$\Delta \mathcal{S}$
plots for the
$De = 0.2$
case shown in figure 7. Initially, the polymers always respond to the local velocity as they are starting from equilibrium or a completely un-stretched configuration. However, due to the larger extension rate, the intensification from the un-stretched configuration at
$H = 0$
to
$0.5$
and then through
$2.0$
shown by figure 7 for
$De = 0.2$
is more rapid for
$De = 0.4$
. This continues to larger strains, represented by
$H = 2.85$
in this case of
$De = 0.4$
, where the intensity of the highly stretched polymers around the extensional axis increases even more. Due to a larger polymer relaxation time, once the polymers are stretched by
$\Delta \textit{De}_{\textit{local}} \gt 0$
regions around the extensional axis, they take a longer time to relax from their highly stretched state. This is particularly evident around the extensional axis as the flow speed in that region is large and increases along the streamline. Therefore, the region of stretched polymers forms an elongated wake that extends to a larger downstream distance from the particle than for
$De = 0.2$
discussed earlier. The memory effects are not particularly pronounced in the
$\Delta \mathcal{S} \gt 0$
region around the compressional axis near the front stagnation point because the flow speed in that region is low and decreases along the streamlines. Hence, the polymers have enough time in this low-speed flow to relax and respond only to the local stretching. However, upon increasing
$H$
from 2.85 to 10, the intensity of the region of extra polymer stretch around the compressional axis increases as observed by comparing figures 8(a) and 8(b). The highly stretched region around the extensional axis starts at roughly the same location along the axis for
$De = 0.2$
and
$De = 0.4$
at all
$H$
because just before this region, there is a stagnation point, and the flow speed is very low. Therefore, the start of extensional stretching is not much affected by polymer relaxation time. For this
$De = 0.4$
, there are significant regions of polymer collapse (
$\Delta \mathcal{S} \lt 0$
) around the front and rear stagnation points and on the particle surface around
$45^\circ$
from the axes. These coincide with the regions where
$\Delta \textit{De}_{\textit{local}} \lt 0$
. They are more larger intense at
$De = 0.4$
than
$De = 0.2$
because the undisturbed polymer stretch is larger in the former.
Figure 8. Value of
$\Delta \mathcal{S}$
for
$De = 0.4$
and
$L = 100$
for Hencky strain,
$H$
: (a) 2.85 and (b) 10.0. Compared with the currently undisturbed polymers, in the red regions the polymers are stretched more, and in the blue regions they are collapsed.
The features of
$\Delta \mathcal{S}$
just observed for
$De=0.2$
and
$0.4$
explain the rheological results of §§ 4.1.1 and 4.1.2. Since the intensity of the highly stretched polymers around the extensional axis is more than the collapsed polymers anywhere at each
$H$
, a positive PIPS or
$\hat{\varPi}^{\textit{PP}}$
, is obtained for
$De=0.2$
. An intense wake of highly stretched polymers explains the positive PIPS in figure 5(a) for
$De=0.4$
. Here PIPS initially increases rapidly with
$H$
because the wake intensifies from when it first appears as a highly stretched region around the extensional axis with a shape similar to plots of figure 7. While still remaining very intense, this wake loses some of its thickness from
$H=2.85$
to 10.0 (the red region around the z axis in figure 8
b is slightly thinner than figure 8
a). Yet PIPS increases monotonically between
$H=2.85$
and 10.0 (figure 5
a), because the highly stretched region around the compressional axis intensifies with
$H$
. Dominance of stretching regions over the collapsed regions at all
$H$
explains positive
$\boldsymbol{\hat{S}}{}^{\textit{Vol}}/\hat{\varPi}{}^{U}=\boldsymbol{\hat{S}}{}^{\textit{PP}}/\hat{\varPi}{}^{U}-1$
from figure 5(d).
Figure 9. Contours of
$\Delta\mathcal{S}$
for
$De=0.6$
at
$L=100$
(top row) and
$L=10$
(bottom row), shown for different
$\log(L)$
-normalized Hencky strains. Panels correspond to (a), (d)
$H/\log(L)=0.5$
; (b), (e)
$1.5$
; and (c), (f)
$8.0$
. Red regions indicate polymers stretched relative to the current undisturbed state, while blue regions indicate collapsed polymers.
In figure 9, we show
$\Delta \mathcal{S}$
for
$De=0.6$
at
$L=10$
and 100 and three different
$H/\log (L)=0.5$
, 1.5 and 8.0. figure 9(a–c) are for
$L=100$
and figure 9(d–f) are for
$L=10$
. The steady state is represented by
$H/\log (L)=8.0$
. First consider
$L=100$
. Up to
$H/\log (L)\approx0.5$
similar mechanisms are at play as those we described above for
$De=0.4$
and
$L=100$
up to its steady state. This also explains the similar initial variation of PIPS,
$\hat{\varPi}^{\textit{PP}}/\hat{\varPi}{}^{U}$
in the case of
$De=0.6$
and
$L=100$
as the variation of PIPS up to the steady state of
$De=0.4$
and
$L=100$
(compare figures 5
a and 5
b). For the
$De=0.6$
case the effects are more pronounced. The wake of stretched polymers around the extensional axis is more intense (compare figures 8
b and 9
a– both their wakes appear similar but the colour scale for the latter is over a larger range) that leads to larger positive increase in the PIPS (compare figures 5
a and initial part of 5
b). But, unlike
$De=0.4$
, this is not the steady state for
$De=0.6$
.
For
$De=0.6$
, the undisturbed polymers keep getting more stretched with
$H$
and undergo coil-stretch transition where they stretch close to their maximum extensibility
$L$
. Hence, the polymers around the particle surface undergo an intense relative collapse compared with the highly stretched, undisturbed polymers. This is observed by a much larger and more intense blue region in figure 9(b) at
$H/\log (L)=1.5$
as compared with
$H/\log (L)=0.5$
in figure 9(a) or at the steady state
$H=8.0$
for De = 0.4 in figure 8(b). The collapsed polymer region extends further downstream from the surface due to memory effects as the polymers need time and downstream distance to recover from their collapsed state. Initially, the collapse arises because the far-field polymers undergo coil-stretch transition before the polymers near the particle surface have time to respond. The collapse persists into the steady state because when the highly stretched polymers from the far field arrive at the front stagnation point, the low-stretching region represented by blue or
$\Delta \textit{De}_{\textit{local}}\lt 0$
around this location shown in figure 6(a) makes them undergo a stretch-to-coil transition. They cannot fully recover from this collapsed state until a large distance downstream of the particle for
$De=0.6$
and L = 100. The region of collapsed polymers ultimately engulfs the wake of highly stretched polymers around the extensional axis as shown in figure 9(c) for
$H/\log (L)=8.0$
. A more detailed discussion about the steady-state collapse and other features of
$\Delta \mathcal{S}$
in the steady state for a wide range of
$\textit{De}$
and
$L$
is made in Sharma & Koch (Reference Sharma and Koch2023b
).
The beginning of the collapse marks the end of the increasing PIPS in figure 5(b). Beyond this point, PIPS starts to decrease, and due to a large region of collapsed polymers around the particle surface in the steady state represented by figure 9(c), it becomes negative at large
$H / \log (L)$
. The volumetric stresslet
$\boldsymbol{\hat{S}}{}^{\textit{Vol}} / \hat{\varPi}{}^{U} = \boldsymbol{\hat{S}}{}^{\textit{PP}} / \hat{\varPi}{}^{U} - 1$
from figure 5(e) is also negative because of this collapse. These large negative values make the total interaction stress (
$\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{Vol}}$
) negative. Hence, the net negative impact of adding particles discussed in § 4.1.1 arises from this large region of collapsed polymers for
$De = 0.6$
and
$L = 100$
.
At
$L = 10$
, the initial intensification of the wake of highly stretched polymers is weaker, and in the steady state or at large
$H / \log (L)$
, the region of collapse is smaller and less intense than in the
$L = 100$
case. This can be observed by comparing figures 9(a) to 9(c) with the corresponding
$H / \log (L)$
from figures 9(d) to 9( f). A less intense wake in the initial phase for
$L = 10$
implies a smaller increase in
$\hat{\varPi}^{\textit{PP}} / \hat{\varPi}{}^{U}$
shown in figure 5(b). A smaller and less intense region of collapsed polymers implies a smaller and less rapid decrease of
$\hat{\varPi}^{\textit{PP}} / \hat{\varPi}{}^{U}$
in figure 5(b) and a less negative
$\boldsymbol{\hat{S}}{}^{\textit{Vol}} / \hat{\varPi}{}^{U} = \boldsymbol{\hat{S}}{}^{\textit{PP}} / \hat{\varPi}{}^{U} - 1$
in figure 5(e). Therefore, the rheological changes shown in figure 4(c) for the total interaction stress at
$L = 10$
are less profound than at
$L = 100$
for all
$H / \log (L)$
.
Finally, we explain the behaviour of rheological changes observed in §§ 4.1.1 and 4.1.2 at large
$\textit{De}$
through the
$\Delta \mathcal{S}$
field of
$De = 5.0$
shown in figure 10. Initially, a similar mechanism as that previously explained for the initial phase for
$De = 0.6$
occurs. For both
$L = 10$
and 100, as
$H$
increases from 0, a wake of highly stretched polymers forms around the extensional axis, intensifying upon initially increasing
$H$
. Along with this, polymers around the particle surface start to collapse, as shown in figures 10(a) and 10(d) for
$L = 100$
and 10, respectively. This is accompanied by a rapid increase in
$\hat{\varPi}^{\textit{PP}} / \hat{\varPi}{}^{U}$
up to
$H / \log (L) \approx 1.0$
, as shown in figure 5(c).
Figure 10. Contours of
$\Delta\mathcal{S}$
for
$De=5$
at
$L=100$
(top row) and
$L=10$
(bottom row), shown for different
$\log(L)$
-normalized Hencky strains. Panels correspond to (a), (d)
$H/\log(L)=0.5$
; (b), (e)
$1.15$
; and (c), (f)
$2.0$
. Red regions indicate polymers stretched relative to the current undisturbed state, while blue regions indicate collapsed polymers.
The undisturbed polymers at
$De = 5.0$
undergo a local coil-stretch transition starting at
$H / \log (L) \approx 1.0$
, as shown in figure 1(e). At larger
$\textit{De}$
, the undisturbed coil-stretch transition is much more rapid and leads to larger polymer stress, as observed by comparing
$De = 0.6$
, 2.0 and 5.0 plots in figure 1. Polymer stretch follows the same trend in
$\textit{De}$
. In a short time frame from
$H / \log (L) \approx 1.0$
to
$H / \log (L) \approx 1.15$
, the undisturbed polymer stretch increases rapidly to very large values (figure 1
e). Hence, the polymers in the wake around the extensional axis that were highly stretched relative to the undisturbed polymers before the coil-stretch transition (figures 10
a and 10
d) are now less stretched (figures 10
b and 10
e) than the even more stretched undisturbed polymers after a rapid transition. This manifests as a rapid reduction in the intensity of the highly stretched polymer wake around the extensional axis and an increased intensity of the polymer collapse around the particle surface.
A sharp reduction in
$\hat{\varPi}^{\textit{PP}} / \hat{\varPi}{}^{U}$
and
$(\boldsymbol{\hat{S}}{}^{\textit{PP}} + \hat{\varPi}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
in figures 5(c) and 4(e), respectively, to negative values reflects these polymer stretch dynamics in the suspension rheology. The region of collapsed polymers has a large spatial extent due to the memory of the polymers that were previously collapsed in the regions close to the front stagnation line before travelling to the far field along extensional axis. Upon further increasing
$H / \log (L)$
to 2.0 in figures 10(c) and 10( f) for
$L = 100$
and 10, respectively, the region of collapse shrinks towards the particle. This occurs because at large
$\textit{De}$
, a large imposed extension rate allows a quicker recovery of polymers from collapse to attain stretch close to the undisturbed polymer stretch within small distances from the particle surface. This shrinking of the collapsed polymer region appears as a negative peak in
$\hat{\varPi}^{\textit{PP}} / \hat{\varPi}{}^{U}$
for
$De = 2.0$
and 5.0, as shown for the latter in figure 5(c).
At large
$H$
and
$\textit{De}$
, as shown for
$H / \log (L) = 2.0$
and
$De = 5.0$
in figures 10(c) and 10( f) for
$L = 100$
and 10, respectively, polymer collapse is only in a thin region close to the particle surface. The local polymer stretch is the same as the undisturbed stretch elsewhere in the fluid, and the collapse in the thin region is very intense, as
$\hat{\varPi}^{\textit{PP}} / \hat{\varPi}{}^{U}$
beyond the aforementioned negative peak is still negative and significant in figure 5(c). Therefore, the negative effect of the addition of particles for
$De = 5.0$
and
$L = 100$
mentioned at the beginning of § 4.1.1 arises from this thin region of collapsed polymers.
Due to the local polymer stretch being the same as the undisturbed stretch everywhere except in a thin region in the steady state,
$\boldsymbol{\hat{S}}{}^{\textit{Vol}} / \hat{\varPi}{}^{U} = \boldsymbol{\hat{S}}{}^{\textit{PP}} / \hat{\varPi}{}^{U} - 1$
is zero for
$De = 5.0$
(figure 5
f) and 2.0 (not shown) at large
$H$
. Before this steady state, the transient behaviour of
$\boldsymbol{\hat{S}}{}^{\textit{Vol}} / \hat{\varPi}{}^{U}$
is different for
$L = 10$
and 100, as shown in figure 5( f). For larger
$L$
, the higher intensity of the region of collapsed polymers around the particle surface compared with the wake of stretched polymers around the extensional axis (figure 10
a) leads to a negative
$\boldsymbol{\hat{S}}{}^{\textit{Vol}} / \hat{\varPi}{}^{U}$
in the transient. However, for smaller
$L$
, a more intense wake and less intense collapsed polymer region (figure 10
d) lead to a sharp transient increase in
$\boldsymbol{\hat{S}}{}^{\textit{Vol}} / \hat{\varPi}{}^{U}$
with
$H$
before it decreases again.
An analogy between the dynamics of polymer stretch shown in this section and that of a line of dye in a Newtonian fluid is offered in Appendix D.
5. Suspension rheology of concentrated polymeric solutions
In the preceding section, we examined the effects of particle–polymer interactions in a liquid with low polymer concentration,
$c$
, where the influence of polymer configuration on the velocity and pressure fields is minimal. The polymer-induced velocity,
$c\boldsymbol{u}{}^{\textit{PI}}$
(defined as
$\boldsymbol{u} - \boldsymbol{u}^{\textit{Stokes}}$
in (A6)), was small enough that
$| c\boldsymbol{u}{}^{\textit{PI}} | \ll | \boldsymbol{u}^{\textit{Stokes}} |$
, allowing the Newtonian velocity field,
$\boldsymbol{u}^{\textit{Stokes}}$
, to predominantly govern the evolution of the polymer configuration. However, as
$c$
increases,
$\mathcal{O}(| c\boldsymbol{u}{}^{\textit{PI}} |) = \mathcal{O}(|\boldsymbol{u}^{\textit{Stokes}}|)$
, and the polymer configuration depends on the total fluid velocity,
$c\boldsymbol{u}{}^{\textit{PI}} + \boldsymbol{u}^{\textit{Stokes}}$
.
To address this complexity, we utilise direct numerical simulations to study the extensional rheology of dilute suspensions of spheres in concentrated polymeric fluids. These simulations are conducted with our in-house numerical solver. For a detailed methodology, we refer to our previous publication Sharma & Koch (Reference Sharma and Koch2023a ), which includes validation of flow around spheres and spheroids in various viscoelastic flows. We summarise our numerical method in Appendix C, including modifications for evaluating and validating the PIPS.
Here, we first illustrate the rheology results in § 5.1 before discussing the underlying mechanisms in §§ 5.2 and 5.3. The data from the simulations with
$c = 10^{-5}$
(equivalent to that from the low
$c$
methodology as shown in Appendix C) represent the low
$c$
results to be compared with numerical data at larger
$c$
. As mentioned earlier in § 1 and 2, upon increasing
$c$
, the Giesekus model becomes more appropriate than the FENE-P model for polymer solutions. Thus, for the largest
$c$
results shown below, only the Giesekus model is used. For intermediate values of
$c$
, we use both models, where comparison between the two provides a measure of the effect of polymer entanglement in concentrated solutions in addition to their coiling in dilute solutions. While making this comparison, similar to the undisturbed stress in § 3, we ensure that parameters are chosen such that
$\alpha = L^{-2}$
.
5.1. Rheological observations
To understand the effect of polymer concentration on the particle–polymer interaction stress and extensional viscosity,
$\mu _{\textit{ext}}$
(defined in (2.11)), we briefly review the effects at low polymer concentration,
$c$
, from § 4 and our previous publication (Sharma & Koch Reference Sharma and Koch2023b
). Regardless of
$c$
, the undisturbed polymer extensional viscosity,
$0.5c\hat{\varPi}{}^{U}$
, varies linearly with
$c$
at all
$\textit{De}$
. In the low
$c$
regime, the interaction stress contribution to the extensional viscosity,
$0.5c\phi (\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}})$
, is also linear in
$c$
. Thus, the particle–polymer interaction stress normalised with the undisturbed polymer stress,
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
, is independent of
$c$
at low
$c$
.
In § 4, we found that for
$De \lt 0.5$
,
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
increases with Hencky strain,
$H$
, or time. This increase also occurs for
$De \gt 0.5$
initially, but at a faster rate, leading to a larger increase in interaction stress. For
$De \lt 0.5$
, a steady state is reached with a positive interaction stress value. However, for
$De \gt 0.5$
, due to the coil-stretch transition of the undisturbed polymers, as
$H$
increases,
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
decreases and becomes negative after the initial increase. At higher
$\textit{De}$
, this reduction is more rapid and results in a more negative steady-state value, as shown by Sharma & Koch (Reference Sharma and Koch2023b
).
With this low
$c$
rheology in mind, we first summarise the key findings at larger
$c$
before discussing them in more detail using the observations from figures 11 to 17. For
$De \gt 0.5$
, increasing
$c$
significantly impacts the negative
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
found at larger Hencky strains,
$H$
. As
$c$
surpasses a certain threshold dependent on
$\textit{De}$
, and the Giesekus parameter
$\alpha$
, or the FENE-P parameter
$L$
, the interaction stress at larger
$H$
becomes increasingly negative, decreasing at a rate of approximately
$c^2$
. This decrease is more pronounced for Giesekus than for FENE-P liquids and is amplified by higher
$\textit{De}$
, larger
$L$
, and lower
$\alpha$
. These changes in
$\textit{De}$
,
$L$
and
$\alpha$
also increase the undisturbed polymer stress (§ 3). Consequently, at sufficient
$H$
, also the parameter regime where undisturbed polymer stress is larger, the interaction of a dilute concentration of spheres with the polymers counteracts the linear increase in undisturbed polymer stress with
$c$
. In other words, adding a small concentration of spheres leads to an even more effective reduction in suspension stress at larger
$c$
,
$\textit{De}$
,
$L$
,
$1 / \alpha$
and
$H$
, where the undisturbed fluid stress is expected to be higher.
Figure 11. Effect of polymer concentration,
$c$
, on
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
. (a) Time or Hencky strain,
$H$
, evolution for Giesekus liquids with
$\alpha = 0.01$
at
$De = 0.4$
, and (b) comparison of the steady-state values for Giesekus fluids with
$\alpha =0.0004$
with FENE-P liquids at
$L = 50$
.
figure 11(a) shows the evolution of the normalised interaction stress for a Giesekus liquid with
$\alpha = 0.01$
at
$De = 0.4$
for different polymer concentrations,
$c$
. For each
$c$
, the evolution of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
is qualitatively similar to the small
$c$
result, showing an increase towards a positive steady state. However, the magnitude of the steady-state value decreases with increasing
$c$
. This effect is more clearly illustrated in figure 11(b), which shows the steady-state value of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
for a Giesekus fluid with
$\alpha = 0.0004$
(orange dashed curve). An identical effect of
$c$
is also observed for a FENE-P liquid with
$L = 50$
(solid black curve). From this and simulations (not shown) at other values of
$\textit{De}$
,
$L$
and
$\alpha$
at different
$c$
, we find that, for
$De \lt 0.5$
, the interaction stress of a FENE-P liquid with a given
$L$
is similar to that of a Giesekus liquid with
$\alpha = L^{-2}$
when
$L \geqslant 50$
. This observation aligns with the undisturbed polymer stress behaviour across these models (§ 3, figure 2
a).
For
$De \gt 0.5$
, the polymer concentration
$c$
has a more profound influence on
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
, particularly affecting the negative interaction stress values at large
$H$
. This is illustrated for
$De = 2.0$
in figure 12 for Giesekus and FENE-P liquids with
$\alpha = 0.0004$
and
$L = 50$
, and for
$De = 5.0$
in figure 13 with
$\alpha = 0.01$
and
$L = 10$
. The qualitative behaviour of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
evolution with
$H$
is similar across different
$c$
values: an initial increase to large positive values followed by a decrease to large negative values. Similar to the steady-state interaction stress for
$De \lt 0.5$
, the maximum value for
$De \gt 0.5$
slightly reduces with increasing
$c$
.
By comparing the panels in figures 12 and 13, it is evident that this effect on maximum interaction stress is more pronounced for FENE-P than for Giesekus liquids. A clearer indication of this decrease with
$c$
is shown for Giesekus liquids at
$De = 1.0$
with different
$\alpha$
values in figure 14(a), and for
$\alpha = 0.01$
at different
$\textit{De}$
values in figure 14(b).
Figure 12. Effect of polymer concentration,
$c$
, on the evolution of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
at
$De = 2.0$
for (a) Giesekus liquids with
$\alpha = 0.0004$
(simulations running), and (b) FENE-P liquids with
$L = 50$
.
Figure 13. Effect of polymer concentration,
$c$
, on the evolution of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
at
$De = 5.0$
for (a) Giesekus liquids with
$\alpha = 0.01$
, and (b) FENE-P liquids with
$L = 10$
.
Figure 14. Effect of polymer concentration,
$c$
, on the maximum
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
in a Giesekus fluid with (a) different
$\alpha$
at
$De = 1.0$
and (b) with
$\alpha = 0.01$
at different
$\textit{De}$
.
Figure 15. Effect of polymer concentration,
$c$
, on steady-state
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
for a FENE-P (solid lines) and Giesekus (dashed lines) fluid for different
$\textit{De}$
and equivalent
$\alpha$
and
$L$
such that
$\alpha = L^{-2}$
for (a)
$L = 50$
, and (b)
$L = 10$
.
Comparing figures 12(a) and 12(b) reveals a larger magnitude of the negative
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
at larger
$H$
for a Giesekus liquid at
$De = 2.0$
and
$\alpha = 0.0004$
than for a FENE-P liquid at the same
$\textit{De}$
and equivalent
$L = 50$
. Similarly, figures 13(a) and 13(b) at
$De = 5.0$
show that increasing
$c$
makes
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
at larger
$H$
more negative for a Giesekus liquid with
$\alpha = 0.01$
than for a FENE-P liquid with
$L = 10$
. A direct comparison of the effect of
$c$
on steady-state
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
between the two polymer constitutive models at different
$\textit{De}$
and
$L$
values is shown in figure 15 (comparing solid with dashed lines). Therefore, the steady state
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
is more sensitive to
$c$
for Giesekus than FENE-P liquids. In the former, the steady state is approached much slower, consistent with the slower evolution of undisturbed polymer stress in Giesekus fluids shown in figures 2(b) and 2(c) at
$De = 2.0$
and 5.0 with
$L = \alpha ^{-0.5}$
. However, the steady-state
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
is more negative for the Giesekus model, even though the steady-state undisturbed stress of the two models are identical. As shown by the curves in figures 15 and 16, the steady-state normalised interaction stress first reduces in magnitude with
$c$
but then becomes increasingly negative with
$c$
. This latter decrease with
$c$
is almost linear for both FENE-P and Giesekus models.
Comparing the two models for polymer-induced drag reduction in wall-bounded flows, Housiadas & Beris (Reference Housiadas and Beris2013) showed that similar drag reduction is achieved for the Giesekus and FENE-P models when
$L = \alpha ^{-0.5}$
. In contrast, under these conditions, we observe a significant difference in the steady-state interaction stress (and hence extensional viscosity) between the two models at large
$\textit{De}$
. Even when the steady state is reached in the suspension, the individual polymers will continue to stretch and relax transiently as they traverse the steady-state velocity streamlines around the sphere and experience different local flows along their trajectories. While the steady-state interaction stresses in the two models differ significantly at a large Deborah numbers such as
$De=2$
and
$5$
(figures 15
a and 15
b), they are similar to one another at
$De=0.4$
(compare figure 11
b–figure 2
a). This suggests that the differences in the results of the Giesekus and FENE-P models at large
$\textit{De}$
are likely due to the differing coil-stretch transition behaviours of the individual polymers in these constitutive models (figures 2
b and 2
c).
The rate of decrease in steady-state
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
at large
$c$
increases with
$1 / \alpha$
for Giesekus liquids and
$L$
for FENE-P liquids. The influence of
$\alpha$
can be observed by comparing different
$\alpha$
curves at
$De = 1.0$
in figure 16(a) or by comparing the dashed
$De = 1.0$
or 2.0 curves between figures 15(a) and 15(b). Similarly, the effect of
$L$
is observed by comparing different
$L$
curves for
$De = 1.0$
or 2.0 represented by solid lines in figures 15(a) and 15(b). The Deborah number,
$\textit{De}$
, also increases the rate of decrease in steady-state
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
at large
$c$
, as found by comparing different solid lines in figures 15(a) for
$L = 50$
and 15(b) for
$L = 10$
, different dashed lines in figures 15(a) for
$\alpha = 0.0004$
and 15(b) for
$\alpha = 0.01$
, or the different curves in figure 16(b) for a Giesekus fluid with
$\alpha = 0.01$
.
Figure 16. Effect of polymer concentration,
$c$
, on the steady state
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
in a Giesekus fluid with (a) different
$\alpha$
at
$De = 1.0$
and (b) with
$\alpha = 0.01$
at different
$\textit{De}$
.
Figure 17. Effect of polymer concentration,
$c$
, on the variability of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
with
$\textit{De}$
for Giesekus (solid lines) fluids with
$\alpha = 0.1$
and FENE-P (dashed lines) with
$L = 10$
: (a) steady state, and (b) maximum. Both figures share the same legend, and the semi-analytical curve drawn corresponds to a value of −0.853 as estimated by Sharma & Koch (Reference Sharma and Koch2023b
) for FENE-P liquids with small
$c$
.
In our previous work (Sharma & Koch Reference Sharma and Koch2023b
), we found that at large
$\textit{De}$
, the steady-state normalised interaction stress for low
$c$
FENE-P liquids is
$-0.853$
, i.e. independent of
$\textit{De}$
. (We have used slightly different normalising factors across the two studies to define the interaction stress contribution to the suspension stress as found by comparing (2.11) with equation 93 of the previous work.) figure 17(a) shows that simulations of FENE-P liquid (dashed red line with filled circles) with
$c = 10^{-5}$
asymptote to this value at large
$\textit{De}$
. At a given
$c$
, the value that the steady-state normalised interaction stress approaches at large
$\textit{De}$
is different in Giesekus liquids. However, the stress is more sensitive to
$\textit{De}$
at a given
$c$
in Giesekus liquids. At larger
$c$
, the steady-state
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
becomes more sensitive to
$\textit{De}$
, and the value of
$\textit{De}$
at which it plateaus is proportional to
$c$
. As observed in figure 17(b), the maximum value of normalised interaction stress during its evolution increases with
$\textit{De}$
. Here, we can observe that this maximum value is more sensitive for FENE-P than Giesekus liquid.
Figure 18. Effect of polymer concentration,
$c$
, on (a)
$\hat{\varPi}^{\textit{PP}} / \hat{\varPi}{}^{U}$
, and (b)
$\boldsymbol{\hat{S}}{}^{\textit{PP}} / \hat{\varPi}{}^{U}$
at
$De = 2.0$
and
$L = 50$
. Both figures share the same legend.
figure 18 shows the decomposition of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
for different
$c$
for a FENE-P liquid with
$De = 2.0$
and
$L = 50$
into its constituents. The total interaction stress
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
for this case was shown earlier in figure 12(b). We can observe that while the stresslet is quantitatively important, the qualitative changes reported above are primarily driven by PIPS,
$\hat{\varPi}^{\textit{PP}}$
. The details of the qualitative contributions of the stresslet and PIPS to
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
at large
$c$
are the same as those at small
$c$
discussed in § 4.1.2. In §§ 5.2 and 5.3, we will explore the mechanisms driving the rheological changes with
$c$
illustrated in this section.
5.2. Mechanisms reducing the particle–polymer interaction stress magnitude at higher polymer concentrations
As shown in figure 11 for
$De=0.4$
, the interaction stress for
$De\lt 0.5$
reduces at all times upon increasing
$c$
. For
$De\gt 0.5$
, the peak interaction stress decreases with
$c$
, as shown in figure 14. While the aforementioned changes occur for all
$c$
, the magnitude of the steady-state interaction stress at
$De\gt 0.5$
decreases with
$c$
but only up to a moderate
$c$
value that depends on
$\textit{De}$
and
$L$
or
$\alpha$
, as shown in figures 15 and 16. In this section, we investigate the mechanism responsible for this reduction in the effectiveness of the particles to alter the suspension viscosity with
$c$
.
figure 19 illustrates the steady-state extra polymer stretch
$\Delta \mathcal{S}$
around a sphere for
$De = 5.0$
for a Giesekus liquid with different
$c$
values. Increasing
$c$
leads to a reduction in the size of the polymer collapse region and an enhancement of the extra polymer stretch in the surrounding areas. These changes in
$\Delta \mathcal{S}$
contribute to the observed decrease in the magnitude of the negative particle–polymer stress in figure 15 as
$c$
increases from very small values (insets show the numerically calculated values approaching an asymptote expected in the small
$c$
limit). The
$\Delta \mathcal{S}$
plots in figure 19 correspond to the green dashed line in figure 15(b). The correlation between
$\Delta \mathcal{S}$
and
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
remains consistent across other values (not shown) of
$\textit{De}$
and
$\alpha$
considered in figures 15 and 16. This correlation is also similar for FENE-P liquids (not shown).
Figure 19. Steady-state
$\Delta \mathcal{S}$
at
$De = 5.0$
for a Giesekus liquid with
$\alpha = 0.01$
at different polymer concentrations
$c$
. Panels correspond to (a)
$c=10^{-5}$
, (b)
$c=0.2$
and (c)
$c=0.5$
.
The reduction in the magnitude of negative particle–polymer interaction stress at high
$\textit{De}$
and
$H$
with increasing
$c$
is due to the diminished intensity of the underlying polymer collapse region, influenced by the enhanced interaction between fluid velocity and polymer stress at larger
$c$
. figure 20 illustrates the steady-state
$\Delta \textit{De}_{\textit{local}}$
for three different values of
$c$
, corresponding to the steady-state
$\Delta \mathcal{S}$
for Giesekus liquids shown in figure 19.
Polymer collapse in the small
$c$
regime occurs when polymers enter the low-speed, negative
$\Delta \textit{De}_{\textit{local}}$
(blue) region around the stagnation point along the compressional axis, as discussed in § 4.2. In this regime, the negative
$\Delta \textit{De}_{\textit{local}}$
persists due to limited alteration of the flow field by the polymers. However, at larger
$c$
, as the polymers collapse before the steady state (not shown) in this region, the local fluid acts as one with a lower viscosity, allowing for a more rapid flow with enhanced local stretching capabilities. Consequently, in this region, increasing
$c$
results in a rise in steady-state
$\Delta \textit{De}_{\textit{local}}$
to positive values for sufficiently high
$c$
, as shown in the two right panels of figure 20. This shift effectively reduces the extent of the polymer collapse region. Moreover, the changes in negative
$\Delta \textit{De}_{\textit{local}}$
exhibit greater sensitivity to
$c$
at larger
$\textit{De}$
,
$L$
, and
$1 / \alpha$
(not shown), leading to a more pronounced dependence of steady-state interaction stress on
$c$
, as illustrated by the figures discussed in § 5.1.
Figure 20. Steady-state
$\Delta \textit{De}_{\textit{local}}$
at
$De = 5.0$
for Giesekus liquid with
$\alpha = 0.01$
and different polymer concentrations
$c$
. Panels correspond to (a)
$c=10^{-5}$
, (b)
$c=0.2$
and (c)
$c=0.5$
.
The positive steady-state particle–polymer interaction stress observed at small
$\textit{De}$
, as well as the increase to large positive transient values before the undisturbed coil-stretch transition at high
$\textit{De}$
, is attributed to the wake of extra polymer stretch along the extensional axis, as discussed in § 4.2. When
$c$
is increased, the intensity of this wake diminishes (not shown), paralleling the reduction in the intensity of the polymer collapse region in the steady state for high
$\textit{De}$
cases discussed above. This highly stretched wake results from the region of
$\Delta \textit{De}_{\textit{local}} \gt 0$
. Once the polymers start stretching, the local fluid behaves as if it has a higher viscosity than the surrounding fluid, inhibiting flow and reducing the local stretching capability. Thus, as
$c$
increases, the steady-state
$\Delta \textit{De}_{\textit{local}} \gt 0$
becomes less positive (not shown). This reduction limits the extent to which the polymers can be locally stretched, thereby explaining why increasing
$c$
leads to a decrease in interaction stress in these scenarios, as illustrated in figures 11 and 14.
5.3.
Mechanism for increasingly negative interaction stress with
$c$
beyond moderate
$c$
As described in § 5.1 and illustrated in figures 15 and 16, increasing
$c$
beyond a moderate value results in a greater magnitude of the negative particle–polymer interaction stress at large times for
$De \gt 0.5$
. This section investigates the physical mechanisms responsible for these changes as
$c$
increases.
The changes in
$\Delta \mathcal{S}$
and
$\Delta \textit{De}_{\textit{local}}$
near the particle surface, as described in § 5.2, persist as
$c$
increases. If these were the dominant mechanisms, we would expect the magnitude of the steady-state negative interaction stress at high
$\textit{De}$
to continue decreasing with
$c$
, eventually becoming positive as the region of polymer collapse at
$c = 10^{-5}$
transitions to one with extra polymer stretch at
$c = 0.5$
, as illustrated in figure 19. However, as shown in the figures in § 5.1, beyond a certain moderate
$c$
(which depends on
$\textit{De}$
and
$\alpha$
or
$L$
), further increases in
$c$
lead to a progressively more negative value of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
. The underlying cause of these changes occurs around the compressional axis, further upstream of the particle than previously considered. These changes are depicted for suspensions in Giesekus liquids with
$De = 5.0$
and
$\alpha = 0.01$
at three different
$c$
values in figure 21. These values bracket the transition (at
$c=0.1$
) along the dashed green line in figure 15(b) from a reduction in magnitude of the interaction stress to an increase with
$c$
.
Figure 21. Steady-state
$\Delta \mathcal{S}$
at
$De = 5.0$
for Giesekus liquid with
$\alpha = 0.01$
and different polymer concentrations
$c$
in a larger region around the particle than figure 19. Panels correspond to (a)
$c=10^{-5}$
, (b)
$c=0.2$
, and (c)
$c=0.5$
.
At small
$c = 10^{-5}$
, the region around the compressional (radial) axis is characterised by a radially decaying region of extra polymer stretch as shown in figure 21(a) for a Giesekus liquid with
$De=5.0$
and
$\alpha =0.01$
. As
$c$
increases, the portion of this region, more than 3–5 particle radii upstream along the compressional axis transforms into one of polymer collapse (figures 21
b and 21
c). The intensity of this collapse increases with
$c$
(compare figure 21
c with 21
b). Although the change in polymer stretch due to this collapse is smaller than the increase in local stretch closer to the particle (figure 19), it affects a larger volume, outweighing the local effects.
Similar trends in
$\Delta \mathcal{S}$
with
$c$
at other
$\textit{De}$
and
$L$
or
$\alpha$
values (not shown) explain the change in slope of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
vs.
$c$
at moderate
$c$
values in the various curves presented in figures 15 and 16. Compared with Giesekus liquids, FENE-P liquids are more resistant to changes in polymer stretch due to the presence of the particle (not shown) at all polymer concentrations. This leads to a larger slope of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
vs.
$c$
in Giesekus than in FENE-P liquids at an equivalent
$\alpha = L^{-2}$
in figure 15.
figure 22 shows
$\Delta \textit{De}_{\textit{local}}$
for the same parameters as the
$\Delta \mathcal{S}$
plots from figures 21(a)–21(c). Comparing plots across and within these figures reveals that the switch in the sign of the extra polymer stretch,
$\Delta \mathcal{S}$
, correlates with the switch in the sign of the extra velocity stretching,
$\Delta \textit{De}_{\textit{local}}$
, around the compressional axis. This correlation is consistent for FENE-P liquids and across other values of
$\textit{De}$
and
$\alpha$
or
$L$
discussed in § 5.1.
Figure 22. Steady-state
$\Delta \textit{De}_{\textit{local}}$
at
$De = 5.0$
for Giesekus liquids with
$\alpha = 0.01$
and different polymer concentrations
$c$
in a larger region around the particle than figure 20. Panels correspond to (a)
$c=10^{-5}$
, (b)
$c=0.2$
, and (c)
$c=0.5$
.
As discussed in § 4.2.1 and shown in figure 6(b), the extra stretching of the polymers at low
$c$
in the region 3–5 particle radii upstream along the compressional axis results from the separation of streamlines created by the dipole disturbance of the Newtonian velocity field. The streamlines of the net flow disturbance, in the steady state, created by the particle in a Giesekus fluid with
$\alpha =0.01$
and different values of
$c$
are shown in figure 23. The leftmost panel, corresponding to
$c=10^{-5}$
, depicts primarily the Newtonian disturbance field and shows the separation of streamlines near the compressional axis. As the value of
$c$
increases to 0.2 and 0.5, the disturbance velocity field shown in the middle and right panels indicates a change in the direction of disturbance velocity. This change in streamlines aligns with the transition in the region near the compressional axis from
$\Delta \textit{De}_{\textit{local}}\gt 0$
at
$c=10^{-5}$
to
$\Delta \textit{De}_{\textit{local}}\lt 0$
at
$c=0.2$
and 0.5, as discussed above using figure 22.
Figure 23. Streamlines of the disturbance velocity field for
$De = 5.0$
for Giesekus liquid with
$\alpha = 0.01$
at
$c=$
(a)
$10^{-5}$
, (b)
$0.2$
and (c)
$0.5$
.
The change in the disturbance streamlines further from the particle can be explained by the polymer force field induced near the particle surface. The reversal of streamline direction near the compressional axis arises from the reversal of the dipole disturbance relative to that in the Newtonian velocity field due to the polymer stress. In a linear flow, such as the uniaxial extensional flow considered here, the leading-order far-field disturbance in an unbounded fluid around a particle arises from a force dipole placed at the particle centre. Based on a multipole expansion, far from the particle, the dipole disturbance is given by
\begin{equation} \boldsymbol{u}_{\textit{dipole}}=\boldsymbol{D}:\boldsymbol{\nabla}\boldsymbol{G}, \end{equation}
where
$\boldsymbol{D}$
is the force dipole tensor, and
$\boldsymbol{G}$
is the Green’s function of the Stokes equation. For a sphere in a Newtonian fluid, the dipole tensor is
$20\pi /3 \langle \boldsymbol{e}\rangle$
. In the presence of polymer force,
$\boldsymbol{f}=\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varPi}$
, we can approximately estimate the net force dipole tensor as
\begin{equation} \boldsymbol{D}=\frac {20\pi}{3}\langle \boldsymbol{e}\rangle +F_d^\varPi \langle \boldsymbol{e}\rangle =\frac {20\pi}{3}\langle \boldsymbol{e}\rangle -\int {\rm d}V \frac {(\widehat {\boldsymbol{fr}+\boldsymbol{r f}})}{2}, \end{equation}
where
$F_d^\varPi \langle \boldsymbol{e}\rangle$
is the force dipole tensor induced by the polymer force,
$\boldsymbol{f}=\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varPi}$
. As shown in the above equation,
$F_d^\varPi \langle \boldsymbol{e}\rangle$
is obtained via a volume integral involving
$\boldsymbol{f}$
in the fluid volume. The time evolution of
$F_d^\varPi$
(normalised with
$c$
) for different values of
$c$
and the radial limit in the volume integration,
$r_{\textit{cut-off}}$
, from 2 to 10 is shown in figure 24 for the flow around a sphere in a Giesekus liquid with
$\alpha =0.01$
and
$De=5$
. Different colours and line styles represent different
$r_{\textit{cut-off}}$
and
$c$
, respectively. We observe that the integral converges with
$r_{\textit{cut-off}}$
for each
$c$
, indicating that near-field polymers constitute the majority of the polymer-induced dipole,
$F_d^\varPi \langle \boldsymbol{e}\rangle$
. As evident from the closeness of the curves for
$c=0.2$
and
$0.5$
for a given
$r_{\textit{cut-off}}$
, we note that as
$c$
increases,
$F_d^\varPi /c$
collapses. This partially explains the
$c^2$
scaling of the interaction stress discussed in § 5.1; as the dipole velocity disturbance scales as
$c$
, leading to the resulting polymer stress scaling as
$c^2$
.
Figure 24. Polymer-induced force dipole,
$F_d^\varPi /c$
, for a Giesekus liquid with
$\alpha = 0.01$
and
$De = 5.0$
. Different colours represent various integrating radii cutoffs,
$r_{\textit{cut-off}}$
, and different line styles represent various
$c$
.
The polymer-induced dipole,
$F_d^\varPi$
, is negative at all
$H$
and
$c$
shown in figure 24 and decreases monotonically with
$H$
. Therefore, if
$c$
and
$H$
are sufficiently large, the magnitude of
$F_d^\varPi \langle \boldsymbol{e}\rangle$
exceeds that of the Newtonian dipole,
${20\pi}/{3}\langle \boldsymbol{e}\rangle$
, resulting in a net dipole tensor,
$\boldsymbol{D}$
(5.2), that is negative. In other words, at sufficiently high
$c$
and
$H$
, the net dipole disturbance created by the interaction of the particle and polymers opposes the Newtonian dipole disturbance.
In § 5.2, we noted that as
$c$
increases, the near-field polymers around the particle transition from a collapsed state to a more stretched state relative to the undisturbed polymers (figure 19). While the polymer stretch just off the particle changes with
$c$
, there remains a layer (which becomes thinner with increasing
$c$
) of collapsed polymers at the particle surface at all
$c$
(figure 19). This thin layer of collapsed polymers ensures that the structure of the polymer stress gradients does not qualitatively change with
$c$
. This is depicted by the extensional component of the steady-state
$\boldsymbol{f}=\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varPi}$
for different
$c$
in figure 25. Therefore,
$F_d^\varPi$
, as estimated by the volume integral in the near-field region around the particle, does not change sign with
$c$
.
Figure 25. Steady-state extensional component of the polymeric force induced on the fluid flow,
$\boldsymbol{f}/c=\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varPi}/c$
, for
$De = 5.0$
for Giesekus liquid with
$\alpha = 0.01$
at
$c=$
(a)
$10^{-5}$
, (b)
$0.2$
and (c)
$0.5$
.
As discussed above, at large
$H$
and
$\textit{De}$
, increasing
$c$
causes the polymer stress to alter the velocity field. At larger
$c$
, this field exhibits reduced stretching capability in a region approximately 3–5 particle radii upstream of the particle around the compressional axis. The highly stretched polymers arriving from the far field in this large upstream region, with lower velocity stretching capability, undergo a mild stretch-to-coil transition, resulting in an increasingly negative particle–polymer interaction stress with
$c$
, as discussed in § 5.1. In the next section, we analyse the structure of this wake of collapsed polymers and use it to explain the
$c^2$
scaling of the interaction stress.
5.3.1. Structure and scaling of the wake of polymer collapse
The mechanism discussed above explains the change in slope of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
vs.
$c$
from positive to negative at a moderate
$c$
, which depends on
$\textit{De}$
and
$L$
or
$\alpha$
. This transition occurs due to the formation of a wake of polymer collapse around the compressional axis, extending upstream by three particle radii. The linear negative slope observed from moderate to large
$c$
, or equivalently, the scaling of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}})$
with
$-c^2$
, is attributed to the additional polymer stretch within the wake (as illustrated in figure 21). figure 26 shows the change in polymer stretch,
$\Delta \mathcal{S}$
, for a point on the compressional axis, located 20 particle radii upstream from the particle centre. Here, we observe that
$\Delta \mathcal{S}$
varies linearly with
$c$
across a range of
$\alpha$
and
$\textit{De}$
in Giesekus liquids. Notably, the magnitude of the slope of this variation increases with
$\textit{De}$
and
$1/\alpha$
, mirroring the variation of the slope of
$(\hat{\varPi}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}}) / \hat{\varPi}{}^{U}$
shown earlier in figure 16.
Figure 26. The variation of normalised
$\Delta \mathcal{S}(z = 0, r = 20; c, De, \alpha )$
, representing the polymer stretch on the compressional axis at a distance of 20 particle radii upstream of the particle, within the region of self-similar polymer collapse, for different
$\textit{De}$
and
$\alpha$
in a Giesekus liquid. One curve for FENE-P liquid is also included in each figure. Normalisation is applied to facilitate visualisation, using the magnitude of
$\Delta \mathcal{S}$
corresponding to the first point on each curve.
Figure 27. The variation of
$w(z, r; c, De, \alpha ) ={\Delta \mathcal{S}(z, r; c, De, \alpha )}/{|\Delta \mathcal{S}(z = 0, r; c, De, \alpha )|}$
with distance along the extensional axis,
$z$
(a–d), and the rescaled variable
$z/r$
(e–h) for different values of
$r$
(locations along the compressional direction). The specific parameter sets are as follows: (a), (e)
$De = 1$
,
$\alpha = 0.0004$
,
$c = 0.1$
; (b), ( f)
$De = 1$
,
$\alpha = 0.01$
,
$c = 5.0$
; (c), (g)
$De = 1$
,
$\alpha = 0.01$
,
$c = 1.0$
; and (d), (h)
$De = 5$
,
$\alpha = 0.01$
,
$c = 1.0$
. All panels share the same legend displayed in panel (e).
Furthermore, this wake behaves in a self-similar manner. The top panel of figure 27 illustrates the function
\begin{equation} w(z, r; c, De, \alpha ) = \frac {\Delta \mathcal{S}(z, r; c, De, \alpha )}{|\Delta \mathcal{S}(z = 0, r; c, De, \alpha )|}, \hspace {0.1in} r \gt r_{w\textit{ake-self-similar}}\approx 5-10, \end{equation}
which represents the variation of
$\Delta \mathcal{S}(z, r)$
(normalised by the magnitude of
$\Delta \mathcal{S}$
at the compressional axis,
$z = 0$
) along
$z$
, starting from different locations on the
$r$
axis for various values of
$\textit{De}$
,
$\alpha$
and
$c$
. These curves correspond to the region of polymer collapse upstream of the particle along the compressional axis. In each of the bottom panels of figure 27, the same
$w(z, r; c, De, \alpha )$
variation as the respective top panels is presented, but with the horizontal axis rescaled to
$z/r$
. This rescaling allows the curves for different
$r$
to collapse, with a slightly broader range of
$z/r$
for larger values of
$c$
,
$\textit{De}$
, and
$1/\alpha$
, demonstrating the self-similarity.
The linear scaling of the reduction in polymer stretch with the polymer concentration,
$c$
in the self-similar wake upstream of the particle around the compressional axis translates into a
$c^2$
scaling of the reduction in polymer stress relative to the undisturbed value. Hence, beyond a moderate
$c$
the particle–polymer interaction stress becomes increasingly negative at a rate proportional to
$c^{2}$
.
6. Comparison with previous experiments
Experiments conducted during the SHERE (Shear History Extensional Rheology Experiment) (Hall et al. Reference Hall, McKinley, Erni, Soulages and Magee2009; Soulages et al. Reference Soulages, McKinley, Hall, Magee, Chamitoff and Fincke2010; McKinley et al. Reference McKinley, Haward, Jaishankar and Soulages2011; Jaishankar et al. Reference Jaishankar, Haward, Hall, Magee and McKinley2012; McKinley & Jaishankar Reference McKinley and Jaishankar2013) campaign aboard the International Space Station measured the extensional viscosity of a polymer solution with a dilute concentration of spheres using a filament stretching rheometer (McKinley & Sridhar Reference McKinley and Sridhar2002). The rheometer initially consists of a liquid bridge between two circular plates. The plates are pulled apart at an exponential rate, generating homogeneous uniaxial extension at the centre of the liquid bridge. In some experiments, a pre-shear is imparted before the extension by rotating one of the plates. The Weissenberg number,
$\textit{Wi} = \lambda \dot {\gamma}$
, characterises the pre-shear, where
$\dot {\gamma}$
is the shear rate and
$\lambda$
is the polymer relaxation time. A force transducer on the fixed plate and video recordings of the liquid bridge’s mid-plane diameter enable the estimation of the extensional viscosity.
The experimental campaign comprises three distinct sets: SHERE, SHERE-R and SHERE-II. The SHERE and SHERE-R fluids consist of 0.025 wt % narrow polydispersity high molecular weight polystyrene (the polymer) in oligomeric styrene oil (a Newtonian solvent), with a polymer concentration of
$c = 0.09$
. In the SHERE-II experiments, a 3.5 % volume fraction of 6
$\unicode{x03BC}$
m diameter poly(methyl methacrylate) microspheres is suspended in the same polymer solution used in SHERE and SHERE-R. Although the lowest
$\textit{De}$
in these experiments exceeds the highest reported from the simulations in §§ 4.1 and 5.1, a qualitative comparison remains valid. Comparing results from SHERE/SHERE-R and SHERE-II provides unique experimental evidence of the impact of particle–polymer interactions on extensional rheology. The kinematic viscosity of the SHERE fluid is
$0.036\ \mathrm{m}^2\, \mathrm{s^{-1}}$
, and the polymer relaxation time is approximately
$3\,\mathrm{s}$
, leading to an elasticity number
$El = {De}/{Re} = {\lambda \nu}/{a^2} \sim \mathcal{O}(10^{10}).$
This indicates that inertial effects are negligible and the fluid velocity and pressure fields are quasi-steady, as assumed in this work and discussed in § 2.
figure 28(a) shows the effect of particles at
$De = 15$
and no pre-shear (
$\textit{Wi} = 0$
) as
$H$
evolves from the SHERE experiments. Initially, the particles increase the extensional viscosity (
$\mu _{\textit{ext}}$
, defined in (2.11)), with a subsequent decreasing effect as strain develops. This behaviour qualitatively resembles the numerically obtained (from the method of Sharma & Koch Reference Sharma and Koch2023a
) evolution of particle–polymer interaction stress at
$De \gt 0.5$
, as illustrated in figure 28(b) for a particle volume fraction
$\phi =0.035$
in a FENE-P liquid with
$De=2.0$
and
$L=50$
. This qualitative similarity, where particles reduce extensional viscosity at large
$H$
, was also observed in SHERE experiments at other
$\textit{De}$
values (as shown in figure 29) as well as in the numerical results for small
$c$
in § 4.1 and a general
$c$
in § 5.1. The effect of particles at each Hencky strain is quantified by the relative change in viscosity, defined as
$\mu _{\textit{part}} / \mu _{\textit{fluid}}$
(2.12). For large
$H$
, this relative change decreases linearly with
$H$
for
$De = 11$
, 15, and 18.4, as shown in figure 29(a). This trend persists even with moderate pre-shear (
$\textit{Wi} \ne 0$
), as seen in figure 29(b).
Figure 28. Comparison of the effect of particles on the extensional viscosity of the suspension,
$\mu _{\textit{ext}}$
(2.11), for a dilute suspension of spheres with a volume fraction of 3.5 % in a viscoelastic liquid from: (a) SHERE/SHERE-R and SHERE-II experiments by Hall et al. (Reference Hall, McKinley, Erni, Soulages and Magee2009), Soulages et al. (Reference Soulages, McKinley, Hall, Magee, Chamitoff and Fincke2010), McKinley et al. (Reference McKinley, Haward, Jaishankar and Soulages2011), Jaishankar et al. (Reference Jaishankar, Haward, Hall, Magee and McKinley2012) and McKinley & Jaishankar (Reference McKinley and Jaishankar2013) at
$De = 15$
and
$c=0.09$
; and (b) our numerical simulations using the FENE-P model with
$De=2.0$
,
$L=50$
and
$c=0.1$
.
Figure 29. Relative effect of the particles,
$\mu _{\textit{part}} / \mu _{\textit{fluid}}$
(2.12), in SHERE experiments at (a) three different
$\textit{De}$
values with no pre-shear (
$\textit{Wi} = 0$
); and (b) at three different
$\textit{De}$
values with varying
$Wi$
. In (b), different colours and line styles correspond to different
$Wi$
values, and different
$\textit{De}$
values are indicated by different symbols.
Due to noise within the SHERE data, a linear mixed-effects model is used to assess the statistical significance of the reduction in strain hardening induced by particles at
$H \geqslant 2.25$
. The mixed-effects model is defined as
\begin{align} \frac {\mu _{\textit{part}}}{\mu _{\textit{fluid}}} &= \beta _0 + \beta _H (H^k - 2.25) + \beta _{De} De^k + \beta _{Wi} Wi^k + \gamma _0^m + \gamma _H^m (H - 2.25)\nonumber\\ &\quad + \psi ^k, H \geqslant 2.25. \end{align}
Here,
$\beta _H$
,
$\beta _{De}$
and
$\beta _{Wi}$
are the physical (fixed) parameters modelling the effects of
$H$
above 2.25 across various
$\textit{De}$
and
$Wi$
, while
$\psi ^k$
,
$\gamma _H^m$
and
$\gamma _0^m$
represent Gaussian-distributed random errors. Using
$k \in [1, 120]$
points across
$m \in [1, 12]$
curves in figures 29(a) and 29(b) leads to
$\beta _0 = -0.84$
,
$\beta _H = -0.25$
,
$\beta _{De} = 0.056$
and
$\beta _{Wi} = 0.034$
with 95 % confidence intervals of
$[-1.3, -0.37]$
,
$[-0.29, -0.21]$
,
$[-0.03, 0.09]$
and
$[-0.01, 0.08]$
, respectively. The p-values associated with
$\beta _H$
,
$\beta _{De}$
and
$\beta _{Wi}$
are
$6 \times 10^{-8}$
, 0.002 and 0.1. Thus,
$\beta _H = -0.25$
provides statistically significant evidence of reduced extensional viscosity upon particle addition at large
$H \gt 2.25$
. This effect is independent of pre-shear (
$Wi$
) and extension rate (
$\textit{De}$
), as
$\beta _{De}$
and
$\beta _{Wi}$
are an order of magnitude smaller than
$\beta _H$
. Therefore, the SHERE experiments offer compelling evidence of negative particle–polymer interaction stress in a dilute particle suspension of viscoelastic liquid undergoing uniaxial extensional flow at high
$\textit{De}$
and
$H$
through the mechanisms discussed in §§ 4.2, 5.2, and 5.3.
7. Conclusions
In this paper, we investigate the impact of particle–polymer interactions on the extensional rheology of a dilute suspension of spheres in a polymeric solution using a combination of a semi-analytical method and direct numerical simulations. The particles and polymers engage in a two-way interaction that generates additional stress, referred to as particle–polymer interaction stress, in the suspension. This stress is in addition to the combined stress from a particle-free viscoelastic fluid and the Newtonian particle-induced suspension stress. The latter is
$5\phi \langle \boldsymbol{e}\rangle$
, where
$\phi$
is the particle volume fraction and
$\langle \boldsymbol{e}\rangle$
is the imposed rate of strain tensor. Particles affect the polymers by modifying the flow field and changing their configuration, while the polymers impact the particles by inducing additional elastic traction on the particle surface and altering the local flow field. As a result, a change in the particle-induced extensional viscosity,
$\mu _{\textit{part}}$
(2.12), beyond the Newtonian value of
$2.5\phi$
, is observed. Assuming a dilute particle concentration, i.e. negligible particle–particle interactions, allows us to use ensemble averaging techniques and the flow field around an isolated particle in an unbounded fluid to characterise the suspension rheology.
When the Deborah number (
$\textit{De}$
), defined as the non-dimensional product of polymer relaxation time and the imposed extension rate, is smaller than 0.5, the interaction stress evolves from zero to a positive steady-state value. For
$De \gt 0.5$
, the suspension exhibits a transient behaviour similar to that of lower
$\textit{De}$
suspensions, but with a more pronounced initial increase in interaction stress, which subsequently peaks and then falls to a negative steady state. The steady-state extensional viscosity of particle-free viscoelastic fluids (
$\mu _{\textit{fluid}}$
, as defined in (2.12)) increases with
$\textit{De}$
at large
$\textit{De}$
. This increase results from significant polymer stretching due to the coil-stretch transition (De Gennes Reference De Gennes1974). However, in the
$De \gt 0.5$
regime at large times,
$\mu _{\textit{part}}$
is negative, and thus the suspension viscosity is found to be lower than that of the particle-free fluid. The fluid viscosity,
$\mu _{\textit{fluid}}$
, increases linearly with polymer concentration (
$c$
) and can become substantial in industrial applications. As
$c$
increases, the intensity of particle–polymer interactions rises, and at large
$\textit{De}$
and times, leads to an increased magnitude of the negative interaction stress that scales as
$c^2$
beyond a moderate polymer concentration. Thus, the particles lead to a greater reduction in viscoelastic stress at higher
$c$
and
$\textit{De}$
.
At very low
$c$
, the velocity field in the fluid remains unaffected by the presence of polymers. Thus, using a regular perturbation in
$c$
and a generalised reciprocal theorem, we obtain the suspension rheology up to
$\mathcal{O}(c)$
by using just the analytically available Newtonian velocity field around the particle. These computations are highly efficient as the equations governing mass and momentum conservation in the fluid need not be solved. In this low
$c$
regime, viscoelastic fluids can be characterised as polymer solutions within a Newtonian solvent, well described by the FENE-P constitutive equation (Anna & McKinley Reference Anna and McKinley2008). Thus, we solve the FENE-P equations using the method of characteristics to determine the leading-order polymer configuration. This approach is similar to the calculations performed by Koch et al. (Reference Koch, Lee and Mustafa2016) and Sharma & Koch (Reference Sharma and Koch2023b
); however, due to the extreme stretching of fluid elements along the streamlines around the extensional axis and the transient nature of polymer evolution, additional care is taken, as detailed in Appendix B.
The increase in steady-state interaction stress at
$De \lt 0.5$
and during transients for larger
$\textit{De}$
arises from a wake of highly stretched polymers (relative to those undisturbed by the sphere) around the extensional axis. This wake intensifies until steady state for small
$\textit{De}$
and peaks for larger
$\textit{De}$
. The polymers are stretched in this wake because, near the extensional axis, the disturbances created by the sphere enhance the stretching capability of the underlying velocity field beyond that of the imposed flow. In large
$\textit{De}$
cases, the peak interaction stress occurs just before the polymers far from the sphere transition from a coiled to a stretched state, a mechanism that leads to a significant increase in polymer stress. While the undisturbed polymers in the far field experience substantial stretching, the region near the stagnation point on the sphere’s compressional axis experiences low stretching due to small velocity gradients and low flow speed. Consequently, the polymers close to the particle around the compressional axis undergo collapse, resulting in negative particle–polymer interaction stress at large
$\textit{De}$
and large Hencky strains,
$H$
. In extensional flow rheology,
$H$
(product of extension rate and time) is used as the non-dimensional time scale and represents the strain accumulated in the sample being tested. As the polymers first collapse, they occupy a large volume around the particle. However, since the stretching capability of the velocity gradients is greater downstream, around the extensional axis, these collapsed polymers recover some of their stretch, leading to a reduction in the collapse region. This recovery is reflected in a negative peak of the interaction stress value before it achieves the steady state. These values are consistent with our previous work (Sharma & Koch Reference Sharma and Koch2023b
), where at large
$\textit{De}$
, the interaction stress was estimated to be −0.853
$\phi$
times the undisturbed polymer stress at steady state.
Increasing the polymer concentration (
$c$
) enhances the interaction between local polymer stretching and the velocity field, significantly altering the fluid velocity due to the stress induced by the underlying polymers. While the FENE-P model is suitable up to moderate
$c$
, the Giesekus model is more appropriate for high-concentration polymer solutions. Therefore, we perform direct numerical simulations using our in-house numerical solver (Sharma & Koch Reference Sharma and Koch2023a
) to analyse the flow of FENE-P and Giesekus liquids around a sphere across a wide range of
$\textit{De}$
,
$c$
and the model parameters
$L$
and
$\alpha$
. This method employs finite difference discretisation of the mass and momentum conservation equations, as well as the polymer constitutive equations, and is briefly described and validated in the context of extensional rheology in Appendix C.
At large
$c$
, the effect of
$\textit{De}$
on suspension rheology resembles that observed in the
$\mathcal{O}(c)$
calculations described above. However, as
$c$
increases from very small values, the magnitude of interaction stress decreases at all
$\textit{De}$
and
$H$
because stretched polymers alter the underlying flow. In regions where polymers collapse, the local fluid stress is lower, allowing for swifter flow than the far field if
$c$
is sufficiently large. Consequently, as
$c$
increases, velocity stretching capability diminishes in areas of highly stretched polymers while being enhanced in regions of collapse. This results in the polymer stretch approaching the undisturbed state in both the wake of highly stretched polymers around the extensional axis at low
$\textit{De}$
and up to the peak interaction stress at large
$\textit{De}$
, as well as in the polymer collapse region around the particle at large
$\textit{De}$
and
$H$
.
In Newtonian flow, the region upstream of approximately three particle radii along the compressional axis exhibits additional velocity stretching capability, primarily driven by dipole disturbance. Although this region occupies a large volume, the difference in stretching capability relative to the imposed flow is minimal. Consequently, when
$c$
is small, polymers in this region experience only a slight extra stretch compared with the undisturbed flow, a phenomenon observed at all
$\textit{De}$
and
$H$
. Even at large values of these parameters, the contribution to rheology from this region is overshadowed by that from the polymer collapse near the particle. As
$c$
increases, the polymers begin to resist the Newtonian dipole-induced flow, ultimately transforming it into a region with lower velocity stretching capability. This transformation is attributed to the change in sign of the net force dipole caused by the polymer force field near the particle surface. Thus, with the dipole velocity disturbance acting in the opposite direction to the Newtonian disturbance at large
$\textit{De}$
and
$H$
, the polymers in this region collapse relative to those in the far field. Although upon increasing
$c$
, the intensity of polymer collapse near the particle surface diminishes and may even turn into extra stretch, as mentioned above, beyond a certain
$c$
(dependent on
$\textit{De}$
and
$L$
or
$\alpha$
), the particle–polymer interaction stress becomes increasingly negative with
$c$
. This
$c^2$
scaling of the negative interaction stress arises because the gentle polymer collapse in this large region around the compressional axis scales as
$c$
. At different radial locations (compressional axis coordinate) within the wake of polymer collapse upstream of about three particle radii, the polymer stretch exhibits self-similarity along the extensional axis coordinate
$z$
under the rescaling
$z/r$
.
Partial qualitative experimental validation of the numerically obtained negative particle–polymer interaction stress at large
$\textit{De}$
and
$H$
is provided by previous experiments conducted by McKinley and co-authors on the International Space Station (Hall et al. Reference Hall, McKinley, Erni, Soulages and Magee2009; Soulages et al. Reference Soulages, McKinley, Hall, Magee, Chamitoff and Fincke2010; McKinley et al. Reference McKinley, Haward, Jaishankar and Soulages2011; Jaishankar et al. Reference Jaishankar, Haward, Hall, Magee and McKinley2012; McKinley & Jaishankar Reference McKinley and Jaishankar2013). These experiments utilised a polymer solution at a single concentration of
$c = 0.09$
. Notably, our simulations indicate that the interaction stress varies much more rapidly with
$c$
than the linear increase of stress observed in particle-free fluids. This finding highlights the need for a more extensive experimental campaign involving particle suspensions in both FENE-P and Giesekus liquids across a broader range of industrially relevant concentrations. Future experiments with larger
$c$
could be conducted on Earth, where the fluid would be more viscous and the sagging of the liquid bridge in the filament stretching extensional rheometer due to gravity would be minimised.
One aspect not addressed here is polymer depletion near non-interacting, impenetrable surfaces, as observed in previous experiments (Ausserré et al. Reference Ausserré, Hervet and Rondelez1985; Omari, Moan & Chauveteau Reference Omari, Moan and Chauveteau1989a
,Reference Omari, Moan and Chauveteau
b
; Ausserre et al. Reference Ausserre, Edwards, Lecourtier, Hervet and Rondelez1991). The thickness of these depletion layers increases with the polymer radius of gyration (
$R_g$
) and decreases with shear rate (Ausserre et al. Reference Ausserre, Edwards, Lecourtier, Hervet and Rondelez1991) and
$c$
(Omari et al. Reference Omari, Moan and Chauveteau1989b
). Depletion around suspended particles depends qualitatively on the ratio
$R_g/a$
as larger polymers wrap around the particle and smaller ones perceive the particle as an infinite wall (Fuchs & Schweizer Reference Fuchs and Schweizer2002; Lekkerkerker, Tuinier & Vis Reference Lekkerkerker, Tuinier and Vis2024). For SHERE fluids with molecular weight
$2.25 \times 10^6\,\mathrm{g\, mol^{-1}}$
,
$R_g \approx 80$
nm (Terao & Mays Reference Terao and Mays2004; Jaishankar et al. Reference Jaishankar, Haward, Hall, Magee and McKinley2012). Thus, for particles of radius
$a = 3\,\unicode{x03BC}$
m used in these experiments,
$R_g/a \ll 1$
– a regime in which polymer chains flatten and align parallel to the surface (Doxastakis et al. Reference Doxastakis, Chen, Guzmán and de Pablo2004). In this limit, the depletion layer is thin and negligible at low
$\textit{De}$
for particle sizes and polymeric liquids commonly used in rheology experiments. However, at high
$\textit{De}$
and low
$c$
, collapsed polymer layers – predicted by our simulations to contribute to negative particle–polymer interaction stresses – become extremely thin, and depletion effects may become comparable and non-negligible. These effects could be incorporated in future simulations by modelling slip boundary conditions at the particle surface or spatial non-uniformity in
$c$
, following Mavrantzas & Beris (Reference Mavrantzas and Beris1992), within either the semi-analytical framework (§ 4) for dilute systems or the DNS approach (§ 5) for concentrated polymeric solutions.
The effects of particles on the extensional rheology of a viscoelastic liquid can be either beneficial or harmful in industrial applications. Reducing extensional strain hardening through particle addition could lower power consumption in industrial processes such as hydraulic fracturing and extrusion moulding, where polymeric fluids experience strong extensional flow (§ 1). Conversely, in polymer processing operations like film blowing, strain hardening is crucial for achieving uniform film thickness (Münstedt Reference Münstedt2018). In this context, while the addition of particles can be beneficial, it must be carefully managed to ensure that the accumulated Hencky strain remains below the threshold where large positive interaction stress transitions to a negative value.
Acknowledgements
The first author (A.S.) is currently an employee of NTESS. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC (NTESS), a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525.
Funding
This work was supported by NASA Grant 80NSSC23K034, NSF Grant 2206851 and the U.S. Department of Energy, Office of Science, Advanced Scientific Computing Research (ASCR) Early Career Research Program. This work used the Bridges-2 supercomputer at the Pittsburgh Supercomputing Center through allocations CHM240066 and PHY210025 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program (Boerner et al. Reference Boerner, Deems, Furlani, Knuth and Towns2023), supported by U.S. National Science Foundation Grant nos. 2138259, 2138286, 2138307, 2137603, and 2138296.
Declaration of Interests
The authors report no conflict of interest.
Appendix A. More details on rheology constituents
The particle–polymer interaction stress, given by
${c\phi}(\boldsymbol{\hat{\varPi}}{}^{\textit{PP}} + \boldsymbol{\hat{S}}{}^{\textit{PP}})$
as defined in § 2 and expressed in (2.10), encompasses the influence of fluid traction (from both the solvent and the polymers) on the particle surface, as well as the alteration of polymer stress due to the presence of the particle. The first component, which represents the particle–polymer interaction stresslet, is denoted as
$c\phi \boldsymbol{\hat{S}}{}^{\textit{PP}}$
. The second component, referred to as the PIPS, is represented by
${c\phi}\boldsymbol{\hat{\varPi}}{}^{\textit{PP}}$
. In this section, we provide a detailed mathematical analysis of these two stress components.
A.1. Particle stresslet
In a dilute particle suspension with well-separated particles, inter-particle interaction is rare. Thus, studying the imposed flow (uniaxial extension in our case) around an isolated particle in an unbounded fluid suffices to characterise the suspension rheology. The knowledge of the constitutive relation within the particle is not necessary to obtain the particle stresslet. Using the divergence theorem, Batchelor (Reference Batchelor1970) showed that it is the following area integral of a tensor product of stress on the particle surface, with position vector
$\boldsymbol{r}_p$
,
\begin{equation} \boldsymbol{\hat{S}}(\boldsymbol{\sigma}) = \int _{\boldsymbol{r} = \boldsymbol{r}_p} \text{d}A \hspace {0.1in} \frac {1}{2} [{\boldsymbol{rn}} \boldsymbol{\cdot}\boldsymbol{\sigma} + {\boldsymbol{n}} \boldsymbol{\cdot}\boldsymbol{\sigma} {\boldsymbol{r}}] - \frac {1}{3} \boldsymbol{\delta} ({\boldsymbol{n}} \boldsymbol{\cdot}\boldsymbol{\sigma} \boldsymbol{\cdot}{\boldsymbol{r}}). \end{equation}
The stresslet, being a linear function of stress, can be expressed as a sum of the solvent,
$\boldsymbol{\hat{S}}(\boldsymbol{\tau})$
, and polymeric stresslet,
$c\boldsymbol{\hat{S}}(\boldsymbol{\varPi})$
. However, the solvent stress is perturbed by the polymers, and thus the non-Newtonian nature of the underlying fluid also influences
$\boldsymbol{\hat{S}}(\boldsymbol{\tau})$
. Denoting the polymer-induced solvent stress as
$c\boldsymbol{\tau}{}^{\textit{PI}}$
, such that
\begin{equation} \boldsymbol{\tau} = \boldsymbol{\tau}^{\textit{Stokes}} + c\boldsymbol{\tau}{}^{\textit{PI}}, \end{equation}
with
$\boldsymbol{\tau}^{\textit{Stokes}}$
being the stress on the particle if the fluid were Newtonian, the total solvent stresslet,
$\boldsymbol{\hat{S}}(\boldsymbol{\tau})$
, can be written as the sum of the Stokes (
$\boldsymbol{\hat{S}}(\boldsymbol{\tau}^{\textit{Stokes}})$
) and polymer-induced solvent (
$c\boldsymbol{\hat{S}}(\boldsymbol{\tau}{}^{\textit{PI}})$
) stresslet. In the case of a sphere,
$\boldsymbol{\hat{S}}(\boldsymbol{\tau}^{\textit{Stokes}}) = 5\langle \boldsymbol{e}\rangle V_p$
, as initially shown by Einstein (Reference Einstein1905). The interaction stresslet,
${c\phi}\boldsymbol{\hat{S}}{}^{\textit{PP}}$
, introduced in (2.10), is related to the stresslets arising from
$c\boldsymbol{\hat{S}}(\boldsymbol{\varPi})$
and
$c\boldsymbol{\tau}{}^{\textit{PI}}$
via
\begin{equation} V_p \boldsymbol{\hat{S}}{}^{\textit{PP}} = V_p {\hat{S}}{}^{\textit{PP}} \langle \boldsymbol{e}\rangle = \boldsymbol{\hat{S}}(\boldsymbol{\tau}{}^{\textit{PI}}) + \boldsymbol{\hat{S}}(\boldsymbol{\varPi}). \end{equation}
It can therefore also be termed as the net non-Newtonian stresslet.
Irrespective of the polymer concentration, in an inertia-less fluid, the mass and momentum equations are linear in fluid pressure and velocity, and the polymer stress can be considered a forcing term. This realisation, along with the fluid stress decomposition of (A2), allows the governing fluid mass and momentum (2.1) to be written as the sum of a Stokes and a polymer-induced problem, as demonstrated in Sharma & Koch (Reference Sharma and Koch2023a
). The Stokes problem is the same as the original problem described by (2.1), but without the polymeric stress,
$\boldsymbol{\varPi}$
. These equations are
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^{\textit{Stokes}} = 0, \hspace {0.2in} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau}^{\textit{Stokes}} = 0, \end{equation}
where
$\boldsymbol{\tau}^{\textit{Stokes}} = -p^{\textit{Stokes}} \boldsymbol{\delta} + 2\boldsymbol{e}^{\textit{Stokes}}$
is the Newtonian solvent stress, and
$p^{\textit{Stokes}}$
,
$\boldsymbol{u}^{\textit{Stokes}}$
and
$\boldsymbol{e}^{\textit{Stokes}} = (\boldsymbol{\nabla}\boldsymbol{u}^{\textit{Stokes}} + (\boldsymbol{\nabla}\boldsymbol{u}^{\textit{Stokes}})^{{T}}) / 2$
are the Stokes pressure, velocity and strain rate tensor fields. On these equations, the velocity boundary conditions of the original or complete problem described by (2.1) are imposed (i.e. replace
$\boldsymbol{u}$
with
$\boldsymbol{u}^{\textit{Stokes}}$
in (2.3)). The polymer-induced (labelled with a superscript PI) problem governs a polymer-induced solvent stress,
$\boldsymbol{\tau}{}^{\textit{PI}} = -p^{\textit{PI}} \boldsymbol{\delta} + 2\boldsymbol{e}{}^{\textit{PI}}$
, that is driven by the polymer stress divergence
\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}{}^{\textit{PI}} = 0, \hspace {0.2in} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau}{}^{\textit{PI}} = -\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\varPi}, \end{equation}
subject to zero boundary conditions for the velocity on the particle surface and in the far field. Here,
\begin{align} \begin{split} &c\boldsymbol{u}{}^{\textit{PI}} = \boldsymbol{u} - \boldsymbol{u}^{\textit{Stokes}}, \hspace {0.2in} cp^{\textit{PI}} = p - p^{\textit{Stokes}}, \hspace {0.2in} 2\boldsymbol{e}{}^{\textit{PI}} = \boldsymbol{\nabla}\boldsymbol{u}{}^{\textit{PI}} + (\boldsymbol{\nabla}\boldsymbol{u}{}^{\textit{PI}})^{{T}}, \end{split} \end{align}
where
$\boldsymbol{u}{}^{\textit{PI}}$
and
$p^{\textit{PI}}$
are the polymer-induced velocity and pressure fields.
Since the stress
$\boldsymbol{\tau}{}^{\textit{PI}}$
originates due to
$\boldsymbol{\varPi}$
, the polymer-induced solvent stresslet,
$\boldsymbol{\hat{S}}(\boldsymbol{\tau}{}^{\textit{PI}})$
, must also originate due to
$\boldsymbol{\varPi}$
. In Sharma & Koch (Reference Sharma and Koch2023b
), using a generalised reciprocal theorem, we obtain an alternative decomposition of the total non-Newtonian stresslet arising from the particle–polymer interaction as a function of just the polymer stress field and its undisturbed value,
\begin{equation} V_p \boldsymbol{\hat{S}}{}^{\textit{PP}} = \boldsymbol{\hat{S}}\big(\boldsymbol{\varPi}^U\big) + \boldsymbol{\hat{S}}_{v\textit{olume}}(\boldsymbol{\varPi}; \boldsymbol{\varPi}^U), \end{equation}
where for a particle suspension with each particle’s volume
$V_p$
,
\begin{equation} \boldsymbol{\hat{S}}(\boldsymbol{\varPi}^U) = V_p \boldsymbol{\varPi}^U, \text{and,}\ \boldsymbol{\hat{S}}_{v\textit{olume}}(\boldsymbol{\varPi}, \boldsymbol{\varPi}^U) = \int _{V_f} \text{d}V \hspace {0.1in} (\boldsymbol{\varPi}^U - \boldsymbol{\varPi}) : \boldsymbol{\nabla}\boldsymbol{v}. \end{equation}
For spheres, the particle shape-dependent 3-tensor,
$\boldsymbol{v}$
, is
\begin{equation} v_{\textit{ijk}} = \left (\frac {5}{2r^5} - \frac {5}{2r^7}\right ) r_i r_j r_k + \frac {1}{2r^5} (r_j \delta _{ik} + r_k \delta _{\textit{ij}}) + \left (\frac {1}{2r^5} - \frac {5}{6r^3}\right ) r_i \delta _{jk}. \end{equation}
Physically,
$v_{\textit{ijk}} \langle e \rangle _{jk}$
is the velocity disturbance generated by the presence of the particle under consideration in any imposed strain
$\langle e \rangle _{jk}$
. It can be obtained by Stokes flow solution around a sphere in different straining flows available in textbooks such as Leal (Reference Leal2007). Alternatively, as mentioned by Koch et al. (Reference Koch, Lee and Mustafa2016),
$v_{\textit{ijk}}$
is the decaying Newtonian fluid ‘velocity’ field that extensionally deforms at the particle surface.
The alternative stresslet decomposition of (A7), valid for all values of polymer concentration, lends further insight into the suspension rheology as we discuss in § 4.1.2. The undisturbed part,
$\boldsymbol{\hat{S}}(\boldsymbol{\varPi}^U)$
, is simply the stresslet on a fictitious particle placed in the far field such that the stress on its surface is the same as the undisturbed value. The change in the polymer stress in the fluid volume between the particle surface and the far field, due to the presence of the particle, leads to the stresslet
$\boldsymbol{\hat{S}}_{v\textit{olume}}(\boldsymbol{\varPi}, \boldsymbol{\varPi}^U)$
. Viewing the fluid velocity and pressure through the decomposition into Stokes and polymer-induced parts allows us to obtain deeper mechanistic insights into particle–polymer interaction as considered in §§ 5.2 and 5.3. Additionally, while solving the fluid momentum equation for obtaining the results discussed in § 5, the Stokes component is known analytically, and only the polymer-induced momentum (A5) is evaluated numerically.
A.2. Particle-induced polymer stress
Similar to the particle stresslet, the ensemble-averaged polymer stress,
$c\langle \boldsymbol{\hat{\varPi}}\rangle$
, can also be expressed as a function of relevant flow variables around an isolated particle in an unbounded fluid. However, this requires extra care, as simply replacing it with the volume average of the polymeric stress in the suspension (Jain et al. Reference Jain, Einarsson and Shaqfeh2019) leads to a logarithmic divergence when characterising the rheology of a suspension of spheres in a dilute polymeric liquid. Numerically, this may manifest as absent or slow convergence of the required stress with the domain size, as we found through personal communication with the authors of Jain et al. (Reference Jain, Einarsson and Shaqfeh2019). This divergence was first brought to attention by Koch & Subramanian (Reference Koch and Subramanian2006) for a second-order fluid and later for an Oldroyd-B fluid by Koch et al. (Reference Koch, Lee and Mustafa2016).
By identifying the source of this logarithmic divergence as the linearised polymeric stress, Koch et al. (Reference Koch, Lee and Mustafa2016) obtained the correct form of
$c\langle \boldsymbol{\hat{\varPi}}\rangle$
in terms of flow around an isolated particle. The linearisation of the polymer stress is done about the undisturbed polymer stress field,
$c\boldsymbol{\varPi}^U$
, and velocity field,
$\langle \boldsymbol{u}\rangle$
with the particle-induced velocity disturbance driving the linear perturbations in the polymer stress. We have recently extended the original formulation of Koch et al. (Reference Koch, Lee and Mustafa2016) for Oldroyd-B fluid, implemented for steady-state shear rheology, to investigate the steady-state extensional rheology of the suspension of spheres in dilute FENE-P liquid (Sharma & Koch Reference Sharma and Koch2023b
). This formulation can be applied even to the transient case and to different polymer constitutive models. We briefly present the required formulation for the transient rheology of suspensions in FENE-P and Giesekus liquids and refer to the studies mentioned above for details.
The polymer stress and polymer configuration are a sum of the undisturbed (by the particle’s presence), linear and nonlinear parts
\begin{equation} {\boldsymbol{\varPi}} = {\boldsymbol{\varPi}}^U + {\boldsymbol{\varPi}}^L + {\boldsymbol{\varPi}}^{\textit{NL}}, \hspace {0.2in} {\boldsymbol{\varLambda}} = {\boldsymbol{\varLambda}}^U + {\boldsymbol{\varLambda}}^L + {\boldsymbol{\varLambda}}^{\textit{NL}}, \end{equation}
and the velocity field is the sum of the imposed (or undisturbed) velocity field and a perturbation about this
\begin{equation} \boldsymbol{u} = \langle \boldsymbol{u}\rangle + \boldsymbol{u}'. \end{equation}
The undisturbed polymer stress,
${\boldsymbol{\varPi}}^U$
, is obtained by using the undisturbed velocity,
$\langle \boldsymbol{u}\rangle$
(2.4), in the polymer constitutive model given by (2.5) and (2.6). The linearised polymer stress,
${\boldsymbol{\varPi}}^L$
, is given by
\begin{equation} \boldsymbol{\varPi}^{L} = \begin{cases} \dfrac {1}{De} \bigg(f^{U} {\boldsymbol{\varLambda}}^{L} + \dfrac {\text{tr}({\boldsymbol{\varLambda}}^{L})(f^{U})^2}{L^2} \boldsymbol{\varLambda}^{U}\bigg), & \text{FENE-P} \\[10pt] \dfrac {1}{De} \boldsymbol{\varLambda}^{L}, & \text{Giesekus}, \end{cases} \end{equation}
where the linearised configuration,
$\boldsymbol{\varLambda}^{L}$
, is governed by,
\begin{eqnarray} \begin{split} \frac {\partial {\boldsymbol{\varLambda}}^{L}}{\partial t} +& \langle \boldsymbol{u}\rangle \boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\varLambda}}^{L} - {\langle \boldsymbol{\nabla}\boldsymbol{u}\rangle}^{{T}} \boldsymbol{\cdot}{\boldsymbol{\varLambda}}^{L} - {\boldsymbol{\varLambda}}^{L} \boldsymbol{\cdot}\langle \boldsymbol{\nabla}\boldsymbol{u}\rangle - \boldsymbol{\nabla}{\boldsymbol{u}'}^{{T}} \boldsymbol{\cdot}\boldsymbol{\varLambda}^{U} - \boldsymbol{\varLambda}^{U} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}' = \\ &-\begin{cases} \dfrac {1}{De} \left (f^{U} {\boldsymbol{\varLambda}}^{L} + \dfrac {\text{tr}({\boldsymbol{\varLambda}}^{L})(f^{U})^2}{L^2} \boldsymbol{\varLambda}^{U}\right ), & \text{FENE-P}\\[10pt] \dfrac {1}{De} ((1 - 2\alpha ) {\boldsymbol{\varLambda}}^{L} - \alpha ({\boldsymbol{\varLambda}}^{L} \boldsymbol{\cdot}{\boldsymbol{\varLambda}}^{U} + {\boldsymbol{\varLambda}}^{U} \boldsymbol{\cdot}{\boldsymbol{\varLambda}}^{L})), & \text{Giesekus} \end{cases} \end{split} \end{eqnarray}
Following the derivation of Koch et al. (Reference Koch, Lee and Mustafa2016) and Sharma & Koch (Reference Sharma and Koch2023b
), one obtains that the ensemble average of
$\boldsymbol{\varLambda}^{L}$
and hence
$\boldsymbol{\varPi}^{L}$
is zero. This procedure is described as follows: after taking the ensemble average, the above equation can be recast as a linear operator acting on
$\langle \boldsymbol{\varLambda}^{L}\rangle$
(the ensemble average of
$\boldsymbol{\varLambda}^{L}$
) with zero forcing. Zero forcing occurs because the ensemble average of
$-\boldsymbol{\nabla}{\boldsymbol{u}'}^{{T}} \boldsymbol{\cdot}\boldsymbol{\varLambda}^{U} - \boldsymbol{\varLambda}^{U} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}'$
is zero. The initial condition is
$\langle \boldsymbol{\varLambda}^{L}\rangle = 0$
because initially, the polymers are unstretched and hence have zero stress everywhere. Therefore,
$\langle \boldsymbol{\varLambda}^{L}\rangle = 0$
and
$\langle \boldsymbol{\varPi}^{L}\rangle = 0$
at all times (leading to
$\langle \boldsymbol{\hat{\varPi}}\rangle =\boldsymbol{\hat{\varPi}}{}^{U}+\langle \boldsymbol{\hat{\varPi}}{}^{\textit{NL}}\rangle$
). The volume average of
$\boldsymbol{\varLambda}^{L}$
and hence
$\boldsymbol{\varPi}^{L}$
, however, diverges logarithmically in the steady state as discussed in previous studies. The nonlinear polymer stress,
$\boldsymbol{\varPi}^{\textit{NL}}$
decays faster than
$r^{-3}$
in the far field. In case of transient flows, such as the start-up extensional flow considered here, the divergence is observed beyond a certain time proportional to
$\textit{De}$
(as well as the model parameters
$L$
and
$1/\alpha$
). Thus, removing the linearised stress before approximating the ensemble average with the volume average enables us to express
$c\langle \boldsymbol{\hat{\varPi}}\rangle$
with a converging expression that can be obtained by the solution of an imposed flow around an isolated particle
\begin{equation} c\langle \boldsymbol{\hat{\varPi}}\rangle = c\boldsymbol{\hat{\varPi}}{}^{U} + c{\phi}\boldsymbol{\hat{\varPi}}{}^{\textit{PP}}. \end{equation}
The first term is the polymer stress in the absence of the particles, and the second term is the ensemble average of the nonlinear polymer stress,
$\boldsymbol{\hat{\varPi}}{}^{\textit{NL}}$
, or the PIPS in the suspension, expressed as
\begin{equation} \boldsymbol{\hat{\varPi}}{}^{\textit{PP}} = \frac {1}{V_p} \int _{V_f + V_p} \text{d}V \hspace {0.1in} \boldsymbol{\hat{\varPi}}{}^{\textit{NL}} = \frac {1}{V_p} \int _{V_f + V_p} \text{d}V \hspace {0.1in} \boldsymbol{\hat{\varPi}} - \boldsymbol{\hat{\varPi}}{}^{L} - \boldsymbol{\hat{\varPi}}{}^{U}. \end{equation}
It arises due to the change in the polymer configuration by the presence of the particles and has contributions from both the particle and the fluid phase. Here,
$\phi$
is the particle volume fraction, and
$V_f$
and
$V_p$
are the fluid and the particle volume.
A.3. Small time estimates of rheology
Through a regular expansion in
$t$
of flow variables in the mass, momentum, and constitutive equations, Jain et al. (Reference Jain, Einarsson and Shaqfeh2019) obtained small time estimates for an Oldroyd-B fluid. These estimates are also valid for FENE-P and Giesekus equations, as initially, the polymers are in equilibrium and the effect of finite
$L (\lt \infty )$
and
$\alpha (\gt 0)$
, respectively, is not exhibited in these models at small times. At the linear order in time,
$t$
, as in Jain et al. (Reference Jain, Einarsson and Shaqfeh2019), one finds the PIPS,
$\hat{\varPi}^{\textit{PP}}$
, to be zero and the total interaction stresslet,
$\hat{S}^{\textit{PP}}$
, to be that due to the sphere in a Newtonian fluid with viscosity
$c t / De$
\begin{equation} \hat{\varPi}^{\textit{PP}} \langle \boldsymbol{e}\rangle = \mathcal{O}(t^2), \hspace {0.2in} \hat{S}^{\textit{PP}} \langle \boldsymbol{e}\rangle = 5 \langle \boldsymbol{e}\rangle \frac {ct}{De} + \mathcal{O}(t^2). \end{equation}
With the definitions
\begin{equation} {{\hat{S}}}^{U} = \hat{\varPi}{}^{U}, V_p {{\hat{S}}}^{\textit{Vol}} \langle \boldsymbol{e}\rangle = \hat{{\boldsymbol{S}}}_{v\textit{olume}}(\boldsymbol{\varPi}, \boldsymbol{\varPi}^U), V_p \hat{{S}}^{\boldsymbol{\tau}{\textit{PI}}} \langle \boldsymbol{e}\rangle = \hat{{\boldsymbol{S}}}(\boldsymbol{\tau}{}^{\textit{PI}}), \textrm{and}\ V_p \hat{{S}}^{\boldsymbol{\varPi}} \langle \boldsymbol{e}\rangle = \hat{{\boldsymbol{S}}}(\boldsymbol{\varPi}), \end{equation}
the components of the interaction stresslet from the two decompositions are
\begin{equation} \boldsymbol{\hat{S}}{}^{U} = \boldsymbol{\hat{S}}^{\boldsymbol{\tau}\textit{PI}} = \frac {2}{De} t + \mathcal{O}(t^2), \hspace {0.2in} \boldsymbol{\hat{S}}{}^{\textit{Vol}} = \boldsymbol{\hat{S}}^{\boldsymbol{\varPi}} = \frac {3}{De} t + \mathcal{O}(t^2). \end{equation}
The small time equivalence between individual components of the two stresslet decompositions is fortuitous, for, in general, the individual components are not identical.
Appendix B. Regular perturbation and semi-analytical methodology for rheology of particles in viscoelastic fluids with small polymer concentration
This section outlines the semi-analytical method employed to derive the results discussed in § 4. The approach integrates a regular perturbation expansion in polymer concentration (
$c$
), the method of characteristics, and a generalised reciprocal theorem.
A regular perturbation of the relevant flow variables in
$c$
is expressed as:
$\boldsymbol{\sigma} = \boldsymbol{\sigma}^0 + c\boldsymbol{\sigma}^1 + \mathcal{O}(c^2)$
,
$\boldsymbol{\tau} = \boldsymbol{\tau}^0 + c\boldsymbol{\tau}^1 + \mathcal{O}(c^2)$
,
${p} = {p}^0 + c{p}^1 + \mathcal{O}(c^2)$
,
$\boldsymbol{u} = \boldsymbol{u}^0 + c\boldsymbol{u}^1 + \mathcal{O}(c^2)$
,
$\boldsymbol{\varLambda} = \boldsymbol{\varLambda}^0 + c\boldsymbol{\varLambda}^1 + \mathcal{O}(c^2)$
,
$\boldsymbol{\varLambda}^L = \boldsymbol{\varLambda}^{0L} + c\boldsymbol{\varLambda}^{1L} + \mathcal{O}(c^2)$
,
$\boldsymbol{\varPi} = \boldsymbol{\varPi}^0 + c\boldsymbol{\varPi}^1 + \mathcal{O}(c^2)$
and
$\boldsymbol{\varPi}^L = \boldsymbol{\varPi}^{0L} + c\boldsymbol{\varPi}^{1L} + \mathcal{O}(c^2)$
. At the leading order in
$c$
, momentum and mass conservation are the Newtonian flow equations, as observed from (2.1)–(2.2) or from (A4)–(A6). Thus, at the leading order, the velocity is that around a force- and torque-free sphere suspended in a uniaxial extensional flow of a Newtonian fluid
\begin{eqnarray} &u_i^0 = \begin{cases} \langle {e_{\textit{ij}}}\rangle r_j + \dfrac {5}{2} \left (\dfrac {1}{r^7} - \dfrac {1}{r^5}\right ) \langle {e_{jk}}\rangle r_j r_k r_i - \dfrac {1}{r^5} \langle {e_{\textit{ij}}}\rangle r_j, & r \geqslant 1, \\[8pt] 0, & r \lt 1, \end{cases} \end{eqnarray}
where
$\langle {e_{\textit{ij}}}\rangle$
is defined in (2.4) (Leal Reference Leal2007). We obtain the transient extensional rheology of a dilute suspension in a dilute polymeric fluid up to
$\mathcal{O}(c)$
by using just the velocity field of (B1). In (2.5) and (2.6),
$\boldsymbol{u}^0$
drives
$\boldsymbol{\varLambda}^0$
and
$\boldsymbol{\varPi}^0$
. The undisturbed polymer stress
$c\boldsymbol{\varPi}^U$
is obtained from (2.5) and (2.6) for the FENE-P model with
$\langle \boldsymbol{\nabla}\boldsymbol{u}\rangle = \langle \boldsymbol{e}\rangle$
and truncates at the first order in
$c$
. We obtain the leading-order linearised polymer stress,
$\boldsymbol{\varLambda}^{0L}$
and
$\boldsymbol{\varPi}^{0L}$
, from
$\boldsymbol{u} = \boldsymbol{u}^0$
and
$\langle \boldsymbol{u}\rangle = \boldsymbol{r} \boldsymbol{\cdot}\langle \boldsymbol{e}\rangle$
by using (A11) to (A13). Using
$\boldsymbol{\varLambda}^0$
and
$\boldsymbol{\varPi}^{0L}$
, we obtain the leading-order components of the extensional viscosity required in (2.11)
\begin{equation} \mu _{\textit{ext}} = 1 + 2.5\phi + 0.5c(\hat{\varPi}{}^{U} + {\phi}(\hat{\varPi}^{0\textit{PP}} + \boldsymbol{\hat{S}}{}^{0\textit{PP}} + \mathcal{O}(c))), \end{equation}
where
\begin{equation} \hat{\varPi}^{0\textit{PP}} \langle \boldsymbol{e}\rangle = \frac {1}{V_p} \int _{V_f + V_p} \text{d}V \hspace {0.1in} \boldsymbol{\hat{\varPi}}{}^0 - \boldsymbol{\hat{\varPi}}{}^{0L} - \boldsymbol{\hat{\varPi}}{}^{U}, \text{ and,}\ \boldsymbol{\hat{S}}{}^{0\textit{PP}} = \boldsymbol{\hat{S}}{}^{U} + \boldsymbol{\hat{S}}{}^{0{\textit{Vol}}} = \boldsymbol{\hat{S}}{}^{0\boldsymbol{\tau}\textit{PI}} + \boldsymbol{\hat{S}}{}^{0\boldsymbol{\varPi}}, \end{equation}
and
\begin{align}& \boldsymbol{\hat{S}}{}^{U} = \hat{\varPi}{}^{U}, V_p \boldsymbol{\hat{S}}{}^{0{\textit{Vol}}} \langle \boldsymbol{e}\rangle = \boldsymbol{\hat{S}}_{v\textit{olume}}(\boldsymbol{\varPi}^0; \boldsymbol{\varPi}^U),\nonumber\\& V_p \boldsymbol{\hat{S}}{}^{\boldsymbol{\tau}0\textit{PI}} \langle \boldsymbol{e}\rangle = \boldsymbol{\hat{S}}(\boldsymbol{\tau}{}^{\textit{PI}}), \text{ and},\ V_p \boldsymbol{\hat{S}}{}^{0\boldsymbol{\varPi}} \langle \boldsymbol{e}\rangle = \boldsymbol{\hat{S}}(\boldsymbol{\varPi}^0). \end{align}
Obtaining
$\boldsymbol{\tau}^{0{\textit{PI}}}$
requires a solution of the
$\mathcal{O}(c)$
momentum and mass conservation equation obtained from (A5). This is a partial differential equation and requires sophisticated numerical discretisation such as finite difference, finite volume, finite element or spectral methods. We can circumvent this numerical calculation for the
$c \ll 1$
regime if only
$\boldsymbol{\hat{S}}{}^{0\boldsymbol{\tau}\textit{PI}}$
is known. The entire stress field
$\boldsymbol{\tau}^{0{\textit{PI}}}$
is not required. This is because from
$\boldsymbol{\hat{\varPi}}{}^0$
and
$\boldsymbol{\hat{\varPi}}^{0U}$
(already known without evaluating
$\boldsymbol{\tau}^{0{\textit{PI}}}$
), we can obtain
$\boldsymbol{\hat{S}}{}^{U}$
,
$\boldsymbol{\hat{S}}{}^{0{\textit{Vol}}}$
and
$\boldsymbol{\hat{S}}{}^{0\boldsymbol{\varPi}}$
, thus evaluating
$\boldsymbol{\hat{S}}{}^{0\boldsymbol{\tau}\textit{PI}}$
(from (B3)).
Similar to Koch et al. (Reference Koch, Lee and Mustafa2016) and Sharma & Koch (Reference Sharma and Koch2023b
), we use the method of characteristics to obtain
$\boldsymbol{\varLambda}^0$
and
$\boldsymbol{\varLambda}^{0L}$
from
$\boldsymbol{u}^0$
and
$\langle \boldsymbol{u}\rangle$
. The leading-order partial differential (2.6) and (A13) are converted into ODEs along the streamlines of
$\boldsymbol{u}^0$
and
$\langle \boldsymbol{u}\rangle$
, respectively. In our simulations, we generate these streamlines or the ‘mesh’ through numerical integration. To capture regions with large gradients in
$\boldsymbol{\varLambda}^0$
and
$\boldsymbol{\varLambda}^{0L}$
, we employ a higher streamline density or ‘mesh’ resolution near the surface of the particle and along the extensional axis. The streamlines are initiated at a distance of 5000 particle radii upstream and extend a distance of 50 000 particle radii downstream from the centre of the particle. Additionally, the streamlines extend up to approximately 500 particle radii at a 45
$^\circ$
angle from the extensional axis, covering a significant region around the sphere.
Previous studies using this technique, such as Koch et al. (Reference Koch, Lee and Mustafa2016) and Sharma & Koch (Reference Sharma and Koch2023b
), were concerned with the steady-state values of
$\boldsymbol{\varLambda}^0$
and
$\boldsymbol{\varLambda}^{0L}$
. Here, we require the transient values, and the following augmentations to the numerical implementation are necessary. Initially,
$\boldsymbol{\varLambda}^0 = \boldsymbol{\delta}$
, and
$\boldsymbol{\varLambda}^{0L} = 0$
everywhere in the domain. This implies polymers to be in a stress-free or equilibrium state, as expected at the beginning of the startup of a uniaxial extensional flow. In the following, we refer to the calculation of
$\boldsymbol{\varLambda}^0$
and streamlines of
$\boldsymbol{u}^0$
, but the procedure is similar for
$\boldsymbol{\varLambda}^{0L}$
and streamlines of
$\langle \boldsymbol{u}\rangle$
. We define the locations,
$s(t_n)$
, along each streamline/characteristic where the configuration at time step
$t_n$
,
$\boldsymbol{\varLambda}^0(s(\boldsymbol{x}, t_n); t_n)$
, is available. We perform integration along a particular streamline for a specified small time step
$\Delta t = t_{n+1} - t_n$
to obtain the configuration at the next time step,
$\boldsymbol{\varLambda}^0(s(\boldsymbol{x}, t_{n+1}); t_{n+1})$
. The convected locations remain on the same streamline as the Stokes velocity field is temporally constant. Therefore, using the information only on the concerned streamline at time step
$t_{n+1}$
, we interpolate
$\boldsymbol{\varLambda}^0(s(\boldsymbol{x}, t_{n+1}); t_{n+1})$
along each streamline to obtain the values at the starting locations,
$\boldsymbol{\varLambda}^0(s(\boldsymbol{x}, 0); t_{n+1}) \leftarrow \boldsymbol{\varLambda}^0(s(\boldsymbol{x}, t_{n+1}); t_{n+1})$
. We specify boundary conditions at the point on each streamline which is farthest upstream. At that location, the configuration at each time is the spatially constant undisturbed configuration,
$\boldsymbol{\varLambda}^{0U}(t_{n+1})$
, consistent with the current time,
$t_{n+1}$
. Since the calculation on each streamline is independent, the process described above to calculate
$\boldsymbol{\varLambda}^0$
and
$\boldsymbol{\varLambda}^{0L}$
is massively parallelisable. This parallelisability and the fact that only a simple ODE integration and interpolation are needed make this process computationally efficient.
Appendix C. Methodology for direct numerical simulations of viscoelastic fluid with an arbitrary polymer concentration
This section describes the numerical method utilised to obtain the results presented in § 5. The method involves direct numerical simulations of the equations governing the fluid’s mass and momentum, as well as the polymer constitutive (2.6), which are discussed in detail in Sharma & Koch (Reference Sharma and Koch2023a ) along with validation in several flow scenarios. Since we know the Newtonian component of the velocity and pressure fields (solutions of (A4)) analytically, we solve the polymer induced mass and momentum (A5) for the fluid. Here, we provide a summary of the method, outline the additional calculations required, and present the validation for the PIPS, as defined in § 2 and elaborated upon in Appendix A.2.
Our numerical solver is written in prolate spheroidal coordinates, where the particle is one of the coordinate surfaces. Hence, the shape of any prolate spheroid is exactly modelled, and the no-slip velocity boundary conditions are imposed here. In the current work, a prolate spheroid with an aspect ratio of 1.001 represents a sphere. We varied the aspect ratios in the range 1.01 and 1.0001 and found the results insensitive to this change. The schematic of the computational domain is shown in the left panel, and a rendering of the discretised grid is displayed in the middle and right panels of figure 30. The computational domain is bounded by a no-slip particle surface on the inside and a far-field outer boundary where a uniaxial extensional flow is imposed as the velocity boundary condition. The polymer conformation tensor,
$\boldsymbol{\varLambda}$
, is governed by a first-order hyperbolic equation (2.6). Therefore, we only need to apply one boundary condition at the locations where the characteristics of this equation (velocity streamlines) enter the computational domain. Thus, at the incoming streamlines on the outer boundary, the undisturbed configuration,
$\boldsymbol{\varLambda}^U$
, is used as the boundary condition. Finite difference spatial discretisation is used to represent the fluid mass and momentum conservation as well as the polymer constitutive equations. Fourth, sixth, and eighth-order accurate schemes are used in the radial, polar, and azimuthal directions, respectively. The Schur complement method (Furuichi, May & Tackley Reference Furuichi, May and Tackley2011) is used to solve the quasi-steady mass and momentum equations induced by the polymers (A5). The polymer constitutive equations are expressed in log-conformation form (Fattal & Kupferman Reference Fattal and Kupferman2004; Hulsen, Fattal & Kupferman Reference Hulsen, Fattal and Kupferman2005). The benefit of evolving the logarithm of
$\boldsymbol{\varLambda}$
(as first illustrated by Fattal & Kupferman (Reference Fattal and Kupferman2004)) is that it preserves the positive definiteness of
$\boldsymbol{\varLambda}$
and prevents numerical instabilities at large
$\textit{De}$
.
Figure 30. (a): A schematic representation of the computational domain (two-dimensional slice). The (b) displays the discretised computational domain in its entirety, while the (c) provides a zoomed-in view of a region near the particle surface (red). A higher mesh density is utilised in proximity to the particle surface and along the extensional axis to enhance accuracy in these critical areas.
Higher-order upwinding central (HOUC) schemes of Nourgaliev & Theofanous (Reference Nourgaliev and Theofanous2007) are used to stabilise the convective derivative in the polymer equations. In the FENE-P model, to ensure that polymer stretch,
$\text{tr}(\boldsymbol{\varLambda})$
, at any location in the flow is upper bounded by the prescribed maximum extensibility,
$L$
, the numerical technique of Richter, Iaccarino & Shaqfeh (Reference Richter, Iaccarino and Shaqfeh2010) is used to evolve a bespoke equation for
$f = {L^2}/({L^2 - \text{tr}(\boldsymbol{\varLambda})})$
that ensures
$f \gt 0$
at all locations and times. A second-order accurate Crank–Nicholson method is used to time integrate the polymer constitutive equations.
C.1. Additions to the numerical method of Sharma & Koch (Reference Sharma and Koch2023a ) to calculate PIPS
The method for the equations governing the fluid’s mass and momentum (2.1) and the polymer constitutive relations ((2.5) and (2.6)) was considered in Sharma & Koch (Reference Sharma and Koch2023a
). The results were validated for several examples considering different particle aspect ratios, viscoelastic fluid models, levels of fluid inertia and imposed flows. The only missing piece required for rheological calculations is the evaluation of the PIPS. As discussed in Appendix A.2, to obtain the ensemble-averaged PIPS using a volume average, the linearised polymer stress,
${\boldsymbol{\varPi}}^L$
, must be removed from the extra polymer stress (relative to the undisturbed). This requires solving (A12) and (A13) governing the evolution of
${\boldsymbol{\varPi}}^L$
.
The hyperbolic equation (A13) governing
$\boldsymbol{\varLambda}^{L}$
is evolved alongside and in a similar manner as the polymer constitutive equation (2.6) described in Sharma & Koch (Reference Sharma and Koch2023a
). The linearised polymer configuration,
$\boldsymbol{\varLambda}^{L}$
, is defined everywhere and hence must be evaluated inside the particle in addition to the fluid region outside the particle. The only boundary conditions for this are in the far field, i.e.
$\boldsymbol{\varLambda}^{L}_{\textit{far}} = \boldsymbol{0}$
. A log-conformation formulation is not required or possible for
$\boldsymbol{\varLambda}^{L}$
as it is not a positive definite tensor. We do not experience any stability issues arising during the evolution of (A13). The computational domain for this equation extends to the region within the particle. As seen from its governing equation,
$\boldsymbol{\varLambda}^{L}$
is driven by the undisturbed velocity field
$\langle \boldsymbol{u}\rangle = \boldsymbol{r} \boldsymbol{\cdot}\langle \boldsymbol{e}\rangle$
and the perturbation of the local velocity from the undisturbed,
$\boldsymbol{u}'$
. The spatial discretisation required in the fluid domain for the
$\boldsymbol{\varLambda}^{L}$
equation is similar to that for
$\boldsymbol{\varLambda}$
. Inside the particle,
$\boldsymbol{u}' = -\langle \boldsymbol{u}\rangle$
, and hence only the convective term in (A13) needs to be discretised. This is done in a similar fashion as the discretisation of the convective term in the fluid domain, i.e. using the HOUC schemes of Nourgaliev & Theofanous (Reference Nourgaliev and Theofanous2007). The additional implementation and computational time for evolving
$\boldsymbol{\varLambda}^{L}$
are minimal. We obtain the linearised stress,
$\boldsymbol{\varPi}^{L}$
, from (A12), and now we have all the ingredients necessary to evaluate the PIPS or
$\boldsymbol{\hat{\varPi}}{}^{\textit{PP}}$
from (A15). The volume integration is performed separately in the fluid and particle domain using Gaussian quadratures to obtain
$\boldsymbol{\hat{\varPi}}{}^{\textit{PP}} = \hat{\varPi}^{\textit{PP}} \langle \boldsymbol{e}\rangle$
defined in (A15). From the discussion of the low
$c$
regime in § 4.1, we find that the sphere’s influence on the polymer configuration field extends to large distances downstream of the particle where the extra stretch of the polymers is confined to regions close to the extensional axis. This region becomes thinner in the far field. Therefore, to better resolve it, we cluster more points around the extensional axis.
Several parameters introduced by Sharma & Koch (Reference Sharma and Koch2023a
) dictate the shape of the computational domain and the mesh density. The radius of the sphere (or the minor radius of a non-spherical prolate particle) is fixed at 1, while the outer boundary radius is denoted as
$r_{\textit{far}}$
. As mentioned earlier, the equations are solved in prolate spheroidal coordinates. Consequently, the structured grid within the computational domain is defined by
$N_r$
,
$N_\theta$
and
$N_\phi$
, which represent the number of mesh points in the spheroidal, hyperboloidal and azimuthal directions around a spheroid. For a sphere, these correspond to the number of mesh points in the radial, polar and azimuthal directions. To account for the rapid changes in polymer configuration near the particle surface and along the extensional axis, a non-uniform grid is employed, clustering more mesh points in these critical regions (as shown in the middle panel of figure 30). This clustering is controlled by parameters
$c_1$
(for the radial direction) and
$c_2$
(for the polar direction), as defined in Sharma & Koch (Reference Sharma and Koch2023a
), with values fixed at −1.0 and −1.5, respectively.
C.2. Validation and grid convergence of PIPS
Figure 12 of Sharma & Koch (Reference Sharma and Koch2023a
) validates the interaction stresslet,
$\boldsymbol{\hat{S}}{}^{\textit{PP}}$
, for a sphere in uniaxial extensional flow. The stresslet, evaluated using direct numerical simulations (DNS), demonstrates strong agreement with results from our semi-analytical methodology in Appendix B and the simulations conducted by Jain et al. (Reference Jain, Einarsson and Shaqfeh2019). In this section, we focus on validating the PIPS,
$\hat{\varPi}^{\textit{PP}}$
, against these two sources.
We start by validating the method using a polymer concentration of
$c = 10^{-5}$
in the DNS. We then increased
$c$
to
$10^{-4}$
and observed similar results. The mesh sizes in these validation simulations varied, ranging from 200 to 350 in the radial direction, from 150 to 451 in the azimuthal direction, and from 21 to 25 in the polar direction. The domain size was set between
$r_{\textit{far}} = 100$
and 300. Convergence was achieved for mesh density, time step, and domain size in each case presented. In figure 31, we observe a good match between the semi-analytical (dashed curves) and DNS (solid curves) results across a wide range of
$\textit{De}$
,
$L$
and
$H$
. The small discrepancies are expected, as PIPS is calculated as a volume integral and we cannot achieve the extremely high resolution and large computational domain size of the semi-analytical methodology in our simulations.
Figure 31. Comparison of the evolution of normalised PIPS,
$\hat{\varPi}^{\textit{PP}}/\hat{\varPi}{}^{U}$
, from the two methodologies described in Appendix B and Appendix C. The results from the direct numerical simulations (DNS, Appendix C) at a polymer concentration of
$c = 10^{-5}$
are shown with coloured solid lines, while the results from the semi-analytical methodology (Appendix B) are represented by dashed black lines.
As indicated by the results in § 3, the polymer stress increases with the Deborah number (
$\textit{De}$
), polymer concentration (
$c$
) and polymer extensibility in FENE-P liquids (
$L$
, or equivalently the inverse of mobility in Giesekus fluids,
$1/\alpha$
). This translates to more rapid spatial variation in flow and polymer fields near the particle surface. Therefore, as these parameters are increased, the DNS require a finer mesh resolution to accurately resolve the polymer and fluid velocity fields. For some of the simulations reported in § 5.1, we increased
$N_r$
to approximately 1500 and
$N_\theta$
to approximately 1000 to ensure stable numerical simulations. Since the flow is axisymmetric, a relatively low value of
$N_\phi = 25$
is sufficient. These equations are numerically stiff, and time step convergence is obtained for each numerically stable result, ensuring mesh size convergence.
Figure 32. (a): Fore–aft symmetric flow in a compressional plane across the sphere (red). (b): A schematic illustrating the fore-aft and axisymmetric computational domain employed for some of the larger values of
$c$
,
$L$
and
$1/\alpha$
in § 5.
Figure 33. Three-dimensional versus axisymmetric calculation of normalised PIPS,
$\hat{\varPi}^{\textit{PP}}/\hat{\varPi}{}^{U}$
, for
$De = 1.0$
,
$c = 0.2$
,
$L = 10$
using
$N_r = 500$
and
$N_\theta = 351$
.
Figure 34. Comparison of normalised PIPS,
$\hat{\varPi}^{\textit{PP}}/\hat{\varPi}{}^{U}$
, from our calculations with that of Jain et al. (Reference Jain, Einarsson and Shaqfeh2019) across four different
$\textit{De}$
and
$c = 0.471$
for (a) FENE-P liquids with
$L = 100$
and (b) Giesekus liquids with
$\alpha = 0.001$
. Both figures share the same legend shown on the left figure.
The particle shape and imposed extensional flow are axisymmetric about the extensional axis. Furthermore, the flow is fore–aft symmetric, as shown by the left panel in figure 32. Therefore, we can make the simulations
$2N_\phi$
times faster by simulating a fore–aft and axisymmetric domain shown on the right panel of figure 32. The numerical method of Sharma & Koch (Reference Sharma and Koch2023a
) is slightly altered by imposing a fore-aft symmetric and an axisymmetric boundary condition on the vertical and horizontal axes, respectively. As shown in figure 33 for a FENE-P liquid with
$De = 1.0$
,
$c = 0.2$
, and
$L = 10$
, the normalised PIPS,
$\hat{\varPi}^{\textit{PP}}/\hat{\varPi}{}^{U}$
, from the axisymmetric method is indistinguishable from the original three-dimensional method.
Comparison of
$\hat{\varPi}^{\textit{PP}}/\hat{\varPi}{}^{U}$
, the PIPS normalised with the undisturbed stress, evaluated from our method described above with that from Jain et al. (Reference Jain, Einarsson and Shaqfeh2019) shown in figure 34 shows a good agreement during the initial part of the transition. However, our simulations predict higher values thereafter for the FENE-P liquids (figure 34
a) and during the intermediate period for Giesekus liquids (figure 34
b). This discrepancy is likely due to the neglect of the linearised polymer stress,
${\boldsymbol{\varPi}}^L$
, by Jain et al. (Reference Jain, Einarsson and Shaqfeh2019) as discussed in the remainder of this section.
Figure 35. (a) Extra polymer stress,
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
, (b) linearised polymer stress,
$\hat{{\varPi}}{}^{L}$
and (c) nonlinear polymer stress,
$\hat{{\varPi}}^{\textit{NL}}$
around a sphere (boundary denoted by dashed green curve at
$r^2 + z^2 = 1$
) in a uniaxial extensional flow of a FENE-P liquid with
$De = 0.4$
,
$L = 10$
and
$c = 0.471$
at the steady state (Hencky strain,
$H = 6$
). All three figures share the same colour scale labelled in (c).
Figure 36. Same parameters and caption as figure 35, but showing a larger region around the particle. The dashed black curve at
$z = z_\infty = 75$
and a solid green curve representing
$z = r^{-2} R_0^3$
(
$R_0 = 40$
) are relevant for the volume integral in PIPS.
As defined in (A15) and discussed in Appendix A.2, PIPS,
$\hat{\varPi}^{\textit{PP}}$
, is the volume integral of
$\hat{{\varPi}}{}^{\textit{NL}}= \hat{{\varPi}} - \hat{{\varPi}}{}^{L} - \hat{{\varPi}}{}^{U}$
normalised with particle volume. Due to the symmetry of uniaxial extensional flow and the sphere’s shape, the PIPS integral has the same tensor symmetry as the imposed rate of strain tensor
$\langle \boldsymbol{e}\rangle$
and thus we focus on the extensional (33) component of the various stress fields below. Here,
$\hat{{\varPi}}$
,
$\hat{{\varPi}}{}^{U}$
,
$\hat{{\varPi}}{}^{L}$
and
$\hat{{\varPi}}^{\textit{NL}}$
refer to the 33 component of the stress tensors
$\boldsymbol{\hat{\varPi}}$
,
$\boldsymbol{\hat{\varPi}}{}^{U}$
,
$\boldsymbol{\hat{\varPi}}{}^{L}$
and
$\boldsymbol{\hat{\varPi}}{}^{\textit{NL}}$
respectively. The extra stress in the fluid domain created by the particle relative to the undisturbed or far field from the particle is
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
, and Jain et al. (Reference Jain, Einarsson and Shaqfeh2019) use the particle volume normalised volume integral of this quantity to represent PIPS. In order to replace the ensemble-averaged suspension stress with a volume average, as previously discussed, the linearised stress,
${\hat{\varPi}}{}^{L}$
(which has zero ensemble average) must be subtracted to ensure the convergence of the volume integral.
The three panels of figure 35 show
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
,
$\hat{{\varPi}}{}^{L}$
, and
$\hat{{\varPi}}{}^{\textit{NL}}= \hat{{\varPi}} - \hat{{\varPi}}{}^{L} - \hat{{\varPi}}{}^{U}$
respectively for a FENE-P liquid with
$L = 100$
,
$De = 0.4$
and
$c = 0.471$
, corresponding to the black curve in figure 34(a). The computation domain used in these simulations extends to
$r_{{out}} = 800$
, and the number of mesh points is
$N_r = 1000$
and
$N_\theta = 451$
. figure 36 shows the same quantities as figure 35 but over a much larger region (up to 100 particle radii in each direction). From these figures, we can notice that most of the far-field variation in the extra polymer stress,
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
, is within
$\hat{{\varPi}}{}^{L}$
. Thus, the nonlinear polymer stress
$\hat{{\varPi}}^{\textit{NL}}$
is concentrated much closer to the particle surface.
The simulations conducted by Jain et al. (Reference Jain, Einarsson and Shaqfeh2019) employed a rectangular box with extents of 40 and 40–80 particle radii in the compressional and extensional directions, respectively. This is perhaps an appropriate size if
$\hat{{\varPi}}^{\textit{NL}}$
is used in the PIPS integral. However, from the relatively unbounded simulation conducted here with a spherical domain of 800 particle radii, we can observe a substantial volumetric effect of
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
from figure 36(a). Thus, at around 40–80 particle radii, the extra polymer stress leads to a significant contribution. In our simulations, increasing the computational domain affects
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
near the boundaries of the smaller domain, but not
$\hat{{\varPi}}^{\textit{NL}}$
.
Now we explicitly consider the effect of numerical parameters, such as domain size and mesh resolution, on the spatial distributions of
$\hat{{\varPi}}$
,
$\hat{{\varPi}}{}^{L}$
and
$\hat{{\varPi}}^{\textit{NL}}$
, and their subsequent effects on the volume integrals. We first concentrate on the effect of the spatial extent of the volume integral of
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
in the simulation conducted with a computational domain size
$r_{\textit{far}} = 800$
corresponding to the figures 35 and 36, i.e. a FENE-P liquid with
$L = 100$
,
$De = 0.4$
and
$c = 0.471$
. The volume integral is performed in a region within
$z \leqslant r^{-2} R_0^3$
and
$z \leqslant z_\infty$
, and the parameters
$R_0$
and
$z_\infty$
are varied from 20 to 40 and 300 to 800, respectively. The curve
$z \leqslant r^{-2} R_0^3$
is a streamline of the undisturbed uniaxial extensional flow. Far from the particle, a polymer incoming into the computational domain at the upstream outer boundary (
$\sqrt {z^2 + r^2} = r_{\textit{far}}$
) approximately travels along this curve. figure 37 shows the evolution of PIPS as defined by the integral of
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
in dashed lines and as defined by the integral of
$\hat{{\varPi}}^{\textit{NL}}$
(A15) in solid lines of different colours. In panel (a),
$R_0$
is fixed to 40 and
$z_\infty$
is varied, while in panel (b)
$z_\infty$
is fixed to 600 and
$R_0$
is varied. A dotted black line representing the result from Jain et al. (Reference Jain, Einarsson and Shaqfeh2019) is also included in each plot. Initially, all the curves coincide as the polymers are primarily affected by the region very close to the particle. However, at larger Hencky strains,
$H$
, the dashed curves for the integral of
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
diverge, while the solid curves for the integral of
$\hat{{\varPi}}^{\textit{NL}}$
remain indistinguishable from one another. As
$z_\infty$
or
$R_0$
is increased, the integral of
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
at large
$H$
seemingly approaches that of
$\hat{{\varPi}}^{\textit{NL}}$
, and it remains different from the calculations of Jain et al. (Reference Jain, Einarsson and Shaqfeh2019).
Figure 37. Effect of
$z_\infty$
and
$R_0$
(defined in figure 36
c) on the evaluation of PIPS defined by the volume integral of nonlinear polymer stress,
$\hat{{\varPi}}^{\textit{NL}}$
, (solid lines) and extra polymer stress,
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
(dashed lines) from our simulations for a FENE-P liquid with
$L = 100$
,
$De = 0.4$
and
$c = 0.471$
conducted for
$r_{\textit{far}} = 800$
. The dotted line is from Jain et al. (Reference Jain, Einarsson and Shaqfeh2019). Panels (a) and (b) represent
$R_0 = 40$
and
$z_\infty = 600$
, respectively, with different curves representing different
$z_\infty$
and
$R_0$
labelled in the legend.
Each of the curves presented in figure 37 is converged with respect to the mesh size (as well as time step). This is indicated by figure 38(a), where each type of curve (solid or dashed) is indistinguishable across different mesh sizes used (different colours). The effect of the size of the computational domain,
$r_{{out}}$
, is shown in figure 38(b). As expected, increasing the computational domain improves the
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
integral (dashed lines) by bringing it closer to the integral of
$\hat{{\varPi}}^{\textit{NL}}$
(solid lines). However, solid curves representing the volume integral of
$\hat{{\varPi}}^{\textit{NL}}$
are unaffected by the change in mesh size (all other solid curves are hidden behind the black curve). Therefore, besides being mathematically more appropriate, removing the linearised stress,
$\boldsymbol{\hat{\varPi}}{}^{L}$
, from the extra polymer stress before evaluating the volume integral for PIPS is numerically beneficial in obtaining a faster convergence with mesh size.
Figure 38. (a) Effect of grid size on the
$\hat{{\varPi}}^{\textit{NL}}$
(solid lines) and extra polymer stress,
$\hat{{\varPi}} - \hat{{\varPi}}{}^{U}$
(dashed lines) integral curves with grid 1 corresponding to
$N_r = 1000$
,
$N_\theta = 451$
, grid 2 to
$N_r = 1000$
,
$N_\theta = 551$
and grid 3 to
$N_r = 1000$
,
$N_\theta = 551$
. (b) Effect of computational domain size
$R_{{out}}$
. The simulations are for a FENE-P liquid with
$L = 100$
,
$De = 0.4$
, and
$c = 0.471$
. In these curves,
$R_0 = 40$
and
$z_\infty = 500$
is used, and a dotted curve representing results from Jain et al. (Reference Jain, Einarsson and Shaqfeh2019) is added.
Appendix D. An analogy: Polymers stretch like lines of dye released in the Newtonian flow at previous times
In our previous work (Sharma & Koch Reference Sharma and Koch2023b
), we noted an interesting analogy between the way the presence of a sphere in a uniaxial extensional flow changes the stretch of a line of dye in a Newtonian fluid and the way it changes the polymer stretch shown here through the
$\Delta \mathcal{S}$
field. To quantify the former, we introduced a finite time stretch field, or FTS, defined as
\begin{equation} \text{FTS}(t; \boldsymbol{x}_0) = \frac {1}{2t} \ln \bigg(\frac {1}{\tilde {\lambda}_1(t; \boldsymbol{x}_0)}\bigg), \end{equation}
where
$\tilde {\lambda}_1(t; \boldsymbol{x}_0)$
is the minimum eigenvalue of the Cauchy–Green tensor
\begin{equation} \tilde {\boldsymbol{C}}_{0}^t = \boldsymbol{\nabla} _{\boldsymbol{x}_0} \boldsymbol{\tilde {x}}(t; \boldsymbol{x}_0)^{{T}} \boldsymbol{\cdot}\boldsymbol{\nabla} _{\boldsymbol{x}_0} \boldsymbol{\tilde {x}}(t; \boldsymbol{x}_0), \end{equation}
of the backward time flow map
\begin{equation} \boldsymbol{\tilde {x}}(t; \boldsymbol{x}_0) = \boldsymbol{x}_0 - \int _{0}^{t} \boldsymbol{u}(\boldsymbol{\tilde {x}}(\tau ), \tau ) \text{d} \tau, \end{equation}
generated by the velocity field
$\boldsymbol{u}$
. We refer to Sharma & Koch (Reference Sharma and Koch2023b
) for a detailed derivation and discussion. In the undisturbed uniaxial extensional flow,
$\text{FTS}(t; \boldsymbol{x}_0) = 1$
for all
$\boldsymbol{x}_0$
and
$t$
. The change in the local
$\text{FTS}(t; \boldsymbol{x}_0)$
at location
$\boldsymbol{x}_0$
from this value for a chosen
$t$
,
\begin{equation} \Delta \text{FTS}(t; \boldsymbol{x}_0) = \text{FTS}(t; \boldsymbol{x}_0) - 1, \end{equation}
provides information about the way a fluid element or non-diffusive line of dye released time
$t$
ago is currently stretched. The regions of
$\Delta \text{FTS}(t; \boldsymbol{x}_0) \gt 0$
are the locations where the dye is more stretched, and the regions of
$\Delta \text{FTS}(t; \boldsymbol{x}_0) \lt 0$
appear where the dye is less stretched in the presence of a sphere compared with the dye undergoing stretching in an undisturbed uniaxial extensional flow. We show the
$\Delta \text{FTS}(t; \boldsymbol{x}_0)$
figures in figure 39 for various
$t$
.
Figure 39. The change in the finite time stretch field,
$\Delta {\text{FTS}}(t; \boldsymbol{x}_0)$
, due to the sphere in extensional flow for various
$t=$
(a)
$0.05$
, (b)
$1.25$
, (c)
$2.5$
, (d)
$6.0$
, (e)
$12.0$
and (f)
$40.0$
.
In § 4.2.2, we noted that at a particular
$\textit{De}$
,
$\Delta \mathcal{S}$
undergoes all the changes with
$H$
during its transient that
$\Delta \mathcal{S}$
at a lower
$\textit{De}$
undergoes until the latter’s steady state. However, the changes at the lower
$\textit{De}$
occur in a shorter time than at the larger
$\textit{De}$
. From these
$\Delta {\text{FTS}}(t; \boldsymbol{x}_0)$
figures, we observe that the effect of the sphere on the stretch of a polymer at large
$\textit{De}$
is the same as its effect on the stretching of a non-diffusive line of dye released in the flow at a certain time,
$t$
ago. The analogy is valid up to a
$t$
in proportion with the polymer’s memory or its
$\textit{De}$
.
For small
$\textit{De}$
such as 0.2 with limited memory, the steady-state
$\Delta \mathcal{S}$
field at
$H = 2.0$
in figure 7(b) is analogous to
$\Delta {\text{FTS}}(t; \boldsymbol{x}_0)$
at
$t = 0.05$
(figure 39
a). At
$De = 0.4$
,
$\Delta \mathcal{S}$
first undergoes the same transition as
$De = 0.2$
does up to its steady state. Then, the
$\Delta \mathcal{S}$
field at
$H = 2.85$
(close to steady state) and 10.0 (at steady state) for
$De = 0.4$
shown in figures 8(a) and 8(b) is similar to the
$\Delta {\text{FTS}}(t; \boldsymbol{x}_0)$
field at
$t = 1.25$
(figure 39
b). At
$De = 0.4$
,
$\Delta \mathcal{S}$
does not change much from
$H = 2.85$
to its steady state.
The same qualitative changes as
$\Delta \mathcal{S}$
at
$De = 0.4$
, but within a smaller range of
$H$
, are observed for
$\Delta \mathcal{S}$
at
$De = 0.6$
and
$L = 100$
. However,
$\Delta \mathcal{S}$
at
$De = 0.6$
evolves further, as represented through figures 9(a), 9(b) and 9(c) at
$H / \log (L) = 0.5$
, 1.5 and 8.0. This is similar to the changes in the stretch of the line of dye in a Newtonian fluid from
$t = 1.25$
to 12.0, as shown by the
$\Delta {\text{FTS}}(t; \boldsymbol{x}_0)$
plots at these times in figure 39.
At a very large
$De = 5.0$
and
$L = 100$
, further qualitative changes occur in the
$\Delta \mathcal{S}$
field in the later part of its evolution from
$H / \log (L) = 1.15$
to 2.0, as shown in figures 10(b) and 10(c). These are similar to the qualitative changes in the
$\Delta {\text{FTS}}(t; \boldsymbol{x}_0)$
plots of figure 39 from
$t = 12$
to 40.
At
$De \gt 0.5$
, lower
$L$
stops the analogy at a smaller
$t$
. For example, at
$De = 0.6$
with
$L = 100$
, the evolution of
$\Delta \mathcal{S}$
up to its steady state is analogous to
$\Delta {\text{FTS}}(t; \boldsymbol{x}_0)$
up to
$t = 12$
, as noted above. However, for
$L = 10$
, the steady state
$\Delta \mathcal{S}$
shown in figure 9( f) is analogous to
$\Delta {\text{FTS}}(t; \boldsymbol{x}_0)$
at
$t = 2.5$
. The loss in memory due to limited
$L$
was also noted in Sharma & Koch (Reference Sharma and Koch2023b
). In that case, the analogy was made with
$\Delta {\text{FTS}}(t; \boldsymbol{x}_0)$
at different
$t$
across the steady state
$\Delta \mathcal{S}$
at different
$\textit{De}$
. Here, we have identified the analogy of
$\Delta {\text{FTS}}(t; \boldsymbol{x}_0)$
at different
$t$
with the transient evolution of
$\Delta \mathcal{S}$
at a fixed
$\textit{De}$
.
Open access
