A charged finitely extensible dumbbell model: Explaining rheology of dilute polyelectrolyte solutions
Dmitry Shogin Per Amund Amundsen






 [Also at ]The National IOR Centre of Norway, University of Stavanger, 4036 Ullandhaug, Norway 
 dmitry.shogin@uis.no
Department of Energy Resources, University of Stavanger, 4036 Ullandhaug, Norway



 [Also at ]The National IOR Centre of Norway, University of Stavanger, 4036 Ullandhaug, Norway 

Department of Mathematics and Physics, University of Stavanger, 4036 Ullandhaug, Norway



             


A robust non-Newtonian fluid model of dilute polyelectrolyte solutions is derived from kinetic theory arguments. Polyelectrolyte molecules are modeled as finitely elongated nonlinear elastic dumbbells, where effective charges (interacting through a simple Coulomb force) are added to the beads in order to model the repulsion between the charged sections of polyelectrolyte chains. It is shown that the relative strength of this repulsion is regulated by the electric-to-elastic energy ratio, E, which is one of the key parameters of the model. In particular, E accounts for the intrinsic rigidity of polyelectrolyte molecules and can be used to explain the impact of solvent salinity on polyelectrolyte rheology. 
With two preaveraging approximations, the constitutive equations of the resulting fluid model are formulated in closed form. Material functions predicted by the model for steady shear flow, steady extensional flow, small-amplitude oscillatory shear flow, and start-up and cessation of steady shear flow are obtained and  investigated using a combination of analytical and numerical methods. In particular, it is shown how these material functions depend on E. The two limiting cases of the model – uncharged dumbbells (E=0) and rigid dumbbells (E→∞) – are  included in the analysis. It is found that despite its simplicity, the  model predicts most of experimentally observed rheological features of polyelectrolyte solutions.


A charged finitely extensible dumbbell model: 
 Explaining rheology of dilute polyelectrolyte solutions
    P. A. Amundsen
    December 30, 2023
=======================================================================================================







§ INTRODUCTION


Polymeric liquids – liquids containing extremely long molecules – are of great importance for modern engineering and biotechnological applications. In contrast to liquids consisting of small molecules, polymeric liquids are non-Newtonian: they do not obey the laws of classical fluid mechanics and often behave strictly opposite to what one could expect from "ordinary" fluids. <cit.>
Many synthetic and biological polymers are polyelectrolytes: when dissolved in polar liquids (commonly water), they form polyions – parts of the molecular chain become electrically charged. The charges are of identical sign, which leads to repulsive electric forces acting between different parts of the chain. Such forces make the molecules more rigid and rod-like and therefore mechanically different from their electrically neutral analogs, as described by de Gennes et al. <cit.> 
The physics and chemistry of polyelectrolytes have been a hot research topic during the last few decades, with much attention paid to experiments, theory, and numerical simulations (for details, see the recent review by Muthukumar<cit.> and references therein). At the same time, most of the theories consider only equilibrium properties of polyelectrolyte soltuions; very little progress has been made in understanding and modeling their rheology, which is quite specific. One of the crucial features of polyelectrolytes is that their rheological properties are strongly sensitive to the kind and amount of ions present in the solvent. When the charged parts of a polyion attract counterions – e.g., from a dissolved salt – the repulsive force acting between these parts decreases due to screening; hence, the polyions become more flexible. <cit.> This, of course, strongly affects the rheological behavior of the solution.
Two major fields where flows of polyelectrolyte solutions are of  importance are biological physics <cit.> and enhanced oil recovery.  <cit.> In both cases, transient flows in complex geometries are involved. Such flows typically have both shear and extensional components, which for non-Newtonian fluids are coupled in a nontrivial way. In addition, the salinity of the solvent may vary, making the flows of polyelectrolytes more complicated compared to those of uncharged polymers. To understand the dynamics of such flows – even qualitatively – it is necessary to use advanced tensor models based on microscopic physics<cit.> to adequately describe the forces governing the fluid motion; these models must take the polyelectrolyte nature of the polymer into account. At the same time, such models need to be  simple enough in order to be practically applicable at industrial scales and satisfy the demands of applied technology.
We should also remark that any fluid model designed to work for complex flows must correctly describe the response of the fluid to simple shear and shearfree flows in both steady and transient regimes and predict the well-known rheological properties of polymer solutions, such as shear-thinning, extensional thickening, and existence of stress overshoots at flow start-up.<cit.>
Dumbbell models <cit.> are a relatively simple class of fluid models for polymer solutions, derived from microscopic kinetic theory. Each polymer molecule is thought of as composed of two identical spherical beads connected by a spring (elastic dumbbell) or by simply a rigid rod (rigid dumbbell). The configuration of each dumbbell is completely specified by the connector vector, Q, pointing from one bead to the other. Then, the properties of the molecules are defined by the connector force law, F_c(Q). In particular, the Warner force law <cit.> has been found physically relevant,

    F_c(Q)=H Q1-(Q/Q_0)^2.

[tQ]Qconnector vectormfirst used in Eq. (<ref>)
[tF]F_cconnector forceNfirst used in Eq. (<ref>)
[sQ2]Qdumbbell extensionmfirst used in Eq. (<ref>)
[sQ3]Q_0dumbbell extension limitmfirst used in Eq. (<ref>)
[sH]HWarner spring coefficientN·m^-1first used in Eq. (<ref>)
This law describes a nonlinear spring of finite extensibility: the extension of the spring cannot exceed Q_0, and the spring is Hookean with stiffness H at small extensions.   Finite extensibility proves to be a crucial property for understanding the flow of polymer solutions. The resulting macroscopic fluid model is known as the FENE (finitely elongated nonlinear elastic) dumbbell model. <cit.> A preaveraging approximation (closure) made by Peterlin<cit.> allows to formulate the constitutive equations in closed form, facilitating both analytical investigation and numerical simulations of the fluid model. This updated version of the model is now known as FENE-P ("P" for Peterlin) and is probably the best dumbbell model to date. It is, however, restricted to electrically neutral or weakly hydrolyzed polymers.
The pioneering attempt of constructing a dumbbell model describing polyelectrolyte solutions is that of King and Eisenberg,<cit.> who considered a Hookean dumbbell model modified by the presence of effective charges, interacting through an electrostatic Coulomb force. Then, Dunlap and Leal<cit.> constructed an analogous fluid model based on FENE dumbbells; however, they adopted conformation-dependent friction. No closed-form constitutive equation was derived, and numerical simulations revealed a hysteretic behavior of viscosity and relative extension with respect to flow strength. Ait-Kadi et al.<cit.> formulated a constitutive equation for the model of Dunlap and Leal<cit.> using a conformation tensor approach. Since then, the focus has been kept on implementing advanced numerical methods and improving the modeling of the electric repulsion between the charged parts of  polyelectrolyte molecules.<cit.> In this work, we are aiming at developing an effective phenomenological fluid model, useful for understanding the rheology of polyelectrolyte solutions and for qualitative description of their behavior in complex flows. Such a model must encapsulate all the crucial features of polyelectrolytes: orientability, nonlinearity, finite extensibility, and variable intrinsic rigidity of the molecules, and at the same time be as simple as possible without suffering from pathologies of early kinetic theory models.  For these reasons, we take the successful FENE-P dumbbell model as our starting point. Although screened Coulombic interactions (in the Debye-Hückel approximation) are a proper way to describe the repulsion between charged parts of polyelectrolyte molecules, this would also lead to substantial mathematical complexity. Therefore, we adopt the concept of effective charges, reverting to a simple Coulomb force, as proposed by Dunlap and Leal. <cit.>
The paper is organized as follows: In Sec. <ref>, the underlying assumptions of the new polymer fluid model are specified and the closed-form constitutive equations are derived. In Secs. <ref>-<ref>, we investigate the resulting fluid model. The constitutive equations are used to obtain the material functions predicted by the model, and the properties of these functions are studied. This is done for several standard steady and transient flows: steady shear flow (Sec. <ref>), steady extensional flow (Sec. <ref>), small-amplitude oscillatory shear flow (Sec. <ref>), and start-up and cessation of steady shear flow (Sec. <ref>). Each of Secs. <ref>-<ref> is divided into analysis and discussion parts. Finally, conclusions are presented in Sec. <ref>.
Throughout this paper, SI units are used. Scalars, such as temperature T and shear rate γ̇, are written with lightface italic font; vectors, such as velocity v, with boldface Latin; while second-order tensors, such as rate-of-strain tensor γ̇ and stress tensor τ, with boldface Greek. For the stress tensor, the sign convention of Bird et al. <cit.> is adopted.



§ THE CONSTITUTIVE EQUATIONS


The FENE dumbbell model is assumed applicable for sufficiently dilute polymer solutions so that the polymer molecules interact with the molecules of the Newtonian solvent much more strongly than with each other. The polymer-solvent interactions are typically described by an isotropic Stokes's law with conformation-independent coefficient ζ. <cit.> It has been found that, macroscopically, the polymer contribution to the stress tensor of the solution depends on three parameters: the "ideal-gas" pressure (nkT, where n is the number concentration of  dumbbells, k is Boltzmann's constant, and T is the thermodynamic temperature), the dimensionless nonlinearity parameter,

    b=HQ^2_0/kT,

and a time constant, which we shall find it most convenient to define by

    λ = λ_Q= ζ Q_0^212kT.

[sb]bnonlinearity parameter-defined by Eq. (<ref>)
[sZyl0]λtime constant of C-FENE-P dumbbellssdefined by Eq. (<ref>)
[sZyf]ζhydrodynamic drag coefficientkg·s^-1first used in Eq. (<ref>)
[sn]nnumber concentration of dumbbellsm^-3first used in Eq. (<ref>)
[sk]kBoltzmann's constantJ·K^-1first used in Eq. (<ref>)
[sT]Tthermodynamic temperatureKfirst used in Eq. (<ref>)
[sZyl0h]λ_Htime constant of FENE dumbbellssdefined by Eq. (<ref>)
[sZyl0q]λ_Qtime constant of rigid dumbbellssdefined by Eq. (<ref>)
More commonly, 

    λ_H =3λ_Qb=ζ4H

is used in the literature. Our alternative choice of the time constant will be discussed below.
In order to describe qualitatively the electric repulsion between the charged parts of the polyelectrolyte chain, we assume the beads to carry identical effective charges q, interacting via an  electrostatic Coulomb force. To describe this repulsive interaction, an extra term is added to the FENE-P connector force so that

    F_c = H Q1-(Q/Q_0)^2-q^24πε_0 εQQ^3,

where ε_0 is the permittivity of vacuum and ε is the relative permittivity of the solvent.
[sQy]qeffective bead chargeCfirst used in Eq. (<ref>)
[sZye]εrelative permittivity of the solvent-first used in Eq. (<ref>)
[sZye0]ε_0permittivity of vacuumF·m^-1first used in Eq. (<ref>)
Similar to other models describing dilute polymer solutions, the stress tensor can be written as a sum of independent solvent and polymer contributions. Since the solvent is assumed Newtonian, only the polymer contribution, τ, is of rheological interest. It can be written in two forms – the Kramers form and the Giesekus form, respectively, <cit.> 

    τ   = -n ⟨QF_c⟩+nkTδ, 
    τ   = 14nζ⟨QQ⟩_(1).

[ut]τpolymer contribution to  stress tensorPafirst used in Eq. (<ref>)
[ud]δunit tensor-first used in Eq. (<ref>)
[xZz]⟨⟩configuration-space average-first used in Eq. (<ref>)
[yZ](1)upper-convected time (Oldroyd) derivatives^-1first used in Eq. (<ref>)
Here, the angular brackets denote the configuration-space average,<cit.> δ is the unit tensor, while the subscript "(1)" stands for the upper-convected time derivative, introduced by Oldroyd<cit.> to express the rate of change in tensor properties of a fluid element in a coordinate system deforming with the fluid.
Substituting the modified connector force [Eq. (<ref>)] into the Kramers expression for the stress tensor [Eq. (<ref>)], one gets

    τ = -nH ⟨QQ1-(Q/Q_0)^2⟩+nq^24πε_0 ε⟨QQQ^3⟩+nkTδ.

The configuration distribution function of the dumbbells, needed to calculate the averages in this equation, cannot be obtained directly. Therefore, in order to arrive at a closed-form constitutive equation, the averages must be eliminated from Eq. (<ref>). This is made possible by preaveraging the first two terms on the right-hand side of the equation, which is implemented as follows:

    ⟨QQ1-(Q/Q_0)^2⟩   ≈⟨QQ⟩1-⟨ Q^2⟩ /Q^2_0, 
    ⟨QQQ^3⟩   ≈⟨QQ⟩⟨ Q^2 ⟩^3/2.

Approximation (<ref>) is Peterlin's closure, introduced when formulating the FENE-P dumbbell polymer model. <cit.> The corresponding approximation for the Coulomb term [Eq. (<ref>)] is established by analogous arguments. We propose the abbreviation C-FENE-P for this extended polymer model, where "C" stands for "charged".
To describe the strength of the electric interactions, it is convenient to introduce the dimensionless ratio,

    E=q^2/(4πε_0 ε Q_0) kT,
 
[sE2]Eelectric-to-elastic energy ratio-defined by Eq. (<ref>)
between the characteristic potential energy of the electric repulsion and the thermal energy scale of the dumbbells. Alternatively, one can write

    E=z^2 l_B/Q_0,

where l_B is the Bjerrum length<cit.> and z=q/e the valence. Larger values of E correspond to increased electrostatic repulsion between the beads, i.e., to "stiffer" dumbbells. In the limit E→ 0, the original (uncharged) FENE-P model is recovered, while the dumbbells become rigid as E →∞.
[slb]l_BBjerrum lengthmfirst used in Eq. (<ref>)
[se1]eelementary electric chargeCfirst used in Eq. (<ref>)
[szx]zvalence-first used in Eq. (<ref>)
We further introduce the mean-square relative dumbbell extension, x, as

    x=⟨ Q^2⟩Q_0^2,

[sx]xmean-square relative dumbbell extension-defined by Eq. (<ref>)
which is a measure of dumbbell stretching: x reaches its minimal (non-zero) value at equilibrium and x → 1 at maximal stretching, i.e., when Q→ Q_0.
Making use of the preaveraging approximations   (<ref>) and (<ref>) and notations (<ref>)-(<ref>), Eq. (<ref>) can be rewritten as

    τ = - nH ⟨QQ⟩1-x+nHEb⟨QQ⟩x^3/2+nkTδ.

Having defined the C-FENE-P Z-factor by

    Z=11-x-Eb1x^3/2,
 
[sZ]ZZ-factor-defined by Eq. (<ref>)
one writes Eq. (<ref>) as

    τ =-nHZ⟨QQ⟩+nkTδ.

This is identical in form to the corresponding equation of the FENE-P dumbbell model. <cit.> One can still proceed by taking the Oldroyd derivatives of both sides of the equation, eliminating ⟨QQ⟩_(1) using the Giesekus form of the stress tensor (<ref>) and then eliminating ⟨QQ⟩ using Eq. (<ref>) once more. The result is

    b3 Z τ + λτ_(1)-λ{τ-nkTδ}D_t ln Z=-nkTλγ̇,
 
[st]ttime variablesfirst used in Eq. (<ref>)
[uc]γ̇rate-of-strain tensors^-1defined by Eq. (<ref>)
[tv]vfluid velocitym·s^-1first used in Eq. (<ref>)
where

    γ̇ = (∇ v)+(∇ v)^T

is the rate-of-strain tensor, v is the velocity field of the fluid, t is the time variable, while D_t stands for the material derivative. This result is identical to the constitutive equation of the FENE-P dumbbells. Thus, the difference between the FENE-P and the C-FENE-P models is exclusively the appearance of E in the Z-factor.
Taking the trace of Eq. (<ref>) and making some simple rearrangements, one arrives at

    Zx = 3b( 1-tr τ3nkT).

Combining this with the definition of the Z-factor [Eq. (<ref>)] leads to the following algebraic equation for x:

    11-x-E/b√(x)=Z_FENE,

[sX]Xinverse of x-defined by Eq. (<ref>)
where 

    Z_FENE=1+3b( 1-tr τ3nkT)

is the Z-factor of the original FENE-P dumbbell model.
In order to facilitate the solution of Eq. (<ref>), we replace x with its inverse, X,

    X=1x.

Then,

    11-x=XX-1=1+1X-1,

and Eq. (<ref>) becomes

    (Z_FENE-1)+Eb√(X)=1X-1.

We further introduce a function ℱ of two arguments, s>0 and α≥ 0, as the unique real solution of the equation

    s+α√(y)=1y-1

[yf]FENE"of FENE-P dumbbells"-first used in Eq. (<ref>)
[xF]ℱspecial function-defined by Eq. (<ref>)
[sZya]αdummy argument-first used in Eq. (<ref>)
[sy]ydummy variable-first used in Eq. (<ref>)
[ss]sdummy argument-first used in Eq. (<ref>)
with respect to y. Some properties of this function will be used in the following. In particular, at fixed s, ℱ is  monotonically decreasing with α; and at fixed α, it is monotonically decreasing as s increases. Furthermore, ℱ(s,α)>1 on its domain, with

    lim_s →∞ℱ(s,α)=lim_α→∞ℱ(s,α)=1,

and finally, ℱ(s,0) = 1+1/s. Then, Eq. (<ref>) can be solved for X,

    X = ℱ(Z_FENE-1,E/b),

while from Eq. (<ref>), one obtains

    Z = (Z_FENE-1)ℱ(Z_FENE-1,E/b).

The constitutive equations are therefore completely  formulated through expressions (<ref>), (<ref>), and  (<ref>). 
One observes that the C-FENE-P dumbbell model contains the following four parameters: (nkT), b, λ, and E. The first three are precisely those of the original FENE-P dumbbell model, while E is specific to C-FENE-P and describes the intrinsic rigidity of polyelectrolyte molecules, with larger values of E corresponding to stiffer molecules. It also accounts for the salt-sensitivity of polyelectrolytes, E and the salt concentration in the solvent being inversely related: higher salinity means lower values of E and vice versa.
Finally, we consider the rigid dumbbell limit, E→∞. It follows from Eq.  (<ref>) that y → 1 when α→∞ for a finite s. Thus, ℱ(Z_FENE-1,∞)=1 for finite Z_FENE. This is in agreement with physical expectations: an infinitely strong electric repulsion will extend the spring to the upper limit so that Q→ Q_0 and x → 1. This leads to a polymer fluid model with the following constitutive equations: 

    Z_RDBτ + λτ_(1)-λ{τ-nkTδ}D_t ln Z_RDB=-nkTλγ̇, 
    
        Z_RDB=1-tr τ3nkT.

[yr]RDB"of rigid dumbbells"-first used in Eq. (<ref>)
We shall refer to this as the rigid dumbbell (RDB) polymer model. Another rigid dumbbell model, derived using physical assumptions different from ours, was introduced and investigated earlier by Bird et al.<cit.> The two models share a lot of similarities, but they are not equivalent. A detailed comparison between them lies beyond the scope of this work and shall be discussed elsewhere. In what follows, "RDB" refers to the model formulated by Eqs. (<ref>) and (<ref>). 
Note that neither H nor b appears explicitly in Eqs.  (<ref>) and (<ref>). Moreover, H, and therefore b, is not defined for the rigid dumbbells. As a result, the commonly adopted microscopic time scale λ_H is not applicable in the RDB limit, but λ_Q, defined by Eq. (<ref>), is independent of H and hence provides a universal microscopic time scale for the FENE-P, C-FENE-P, and RDB models. This justifies our choice λ=λ_Q.
In Secs. <ref>-<ref>, we shall investigate and analyze the material functions predicted by the C-FENE-P dumbbell and RDB fluid models for some standard flow regimes. The contribution of the Newtonian solvent to the material functions is well understood; therefore, only the polymer contribution to the material functions will be discussed.


§ STEADY SHEAR FLOW




 §.§ Analysis

Steady shear flow can be locally described by a fluid velocity field given by

    v = [ v_1(x_2)        0        0 ],

at any position x, the stress tensor taking the form

    τ=[ τ_11 τ_12    0; τ_12 τ_22    0;    0    0 τ_33 ],

with 

    τ_(1)=
    -[ 2τ_12  τ_22     0;  τ_22     0     0;     0     0     0 ]γ̇_ 12.

[tx]xposition vectormfirst used in Eq. (<ref>)
[sN1]N_1first normal stress differencePadefined by Eq. (<ref>)
[sN1]N_2second normal stress differencePadefined by Eq. (<ref>)
[sZyV1]Ψ_1first normal stress coefficientPa·s^2defined by Eq. (<ref>)
[sZyV91]Ψ_2second normal stress coefficientPa ·s^2defined by Eq. (<ref>)
[sZyg]γ̇shear rates^-1first used in Eq. (<ref>)
[sZyha]ηnon-Newtonian viscosityPa·sdefined by Eq. (<ref>)
The rate-of-strain tensor has only one independent non-zero component γ̇≡γ̇_ 12=γ̇_ 21, and the three standard steady shear flow material functions – non-Newtonian viscosity (η), first normal stress coefficient (Ψ_1), and second normal stress coefficient (Ψ_2) – are defined, respectively, by

    τ_12    = -η(γ̇)γ̇, 
    
    N_1=τ_11-τ_22   = -Ψ_1(γ̇)γ̇^ 2, 
    
    N_2=τ_22-τ_33   = -Ψ_2(γ̇)γ̇^ 2,

where N_1 and N_2 are the normal stress differences.
Substituting Eqs. (<ref>) and (<ref>) into the constitutive equation (<ref>) yields

    b3Zτ_11 = 2 λτ_12γ̇, 
       b3Z τ_12 = - nkT λγ̇, 
       τ_22 = τ_33 = 0.

It follows from Eqs. (<ref>) and (<ref>) that Ψ_1 is directly proportional to the square of the viscosity, the coefficient of proportionality being independent of E; as seen from Eq.  (<ref>), Ψ_2 vanishes identically,

    Ψ_1(γ̇)    = 2nkTη^2(γ̇), 
    Ψ_2(γ̇)    = 0.

Having eliminated τ_11 from Eqs. (<ref>) and (<ref>), one arrives at the following nonlinear algebraic relation between the shear stress and the shear rate:

    (1+2𝕋_12^ 23) ℱ(3+2𝕋_12^ 2b,Eb) 𝕋_12 = -Λ,

[sTij]𝕋_ijdimensionless stress tensor components-defined by Eq. (<ref>)
where 𝕋_12=τ_12/(nkT) is the dimensionless shear stress and Λ = λγ̇ is the dimensionless shear rate. Equation (<ref>) can be solved numerically to calculate the steady shear flow properties of the C-FENE-P dumbbells for arbitrary values of b and E. 
The influence of E on the relative extension of the dumbbells in steady shear flow is shown in Fig. <ref>(a). In general, a higher value of E leads to larger spring extensions, as expected. This effect is more pronounced at equilibrium and at low-to-medium shear rates. At higher shear rates, the dumbbells are already stretched almost to the upper limit by the flow so that the influence of E becomes small.

The predicted impact of salinity on the non-Newtonian viscosity is shown in Fig.  <ref>(b). The C-FENE-P dumbbells are shear-thinning, which is typical for polymer solutions. The impact of E on viscosity is twofold. First, increasing solvent salinity (decreasing E) at any fixed shear rate leads to a reduction in viscosity. This reduction is largest at small shear rates and decreases as the shear rate increases, vanishing as γ̇→∞. Second, the onset of shear-thinning is shifted towards higher shear rates.
The asymptotic behavior of the viscosity curves governed by Eq. (<ref>) can be studied analytically. At very low shear rates, the viscosity approaches its zero-shear-rate value, 

    η_0 = nkTλℱ(3/b,E/b),

which is highly sensitive to E. In the FENE-P limit (E=0), this result simplifies to

    η_0,FENE=3b+3nkTλ,

and in the RDB limit (E →∞), the zero-shear-rate viscosity is

    η_0,RDB=nkT λ.

[sZyhb]η_0zero-shear-rate viscosityPa ·sfirst used in Eq. (<ref>)
At very high shear rates, the asymptotic behavior of the viscosity curves of the RDB, C-FENE-P (independent of E), and FENE-P models is the same,

    η≈√(32)nkTλΛ^-2/3.

Finally, in the FENE-P and RDB limits, the exact analytical solutions of Eq. (<ref>) can be obtained. The exact analytical solution for the FENE-P dumbbells is given, e.g., by Shogin et al.<cit.> In the RDB limit, Eq. (<ref>) reduces to the following cubic equation for 𝕋_12:

    (1+2𝕋_12^ 23) 𝕋_12 = -Λ.

This equation has one real solution, which can be obtained, e.g., using Cardano's method, as described by La Nave and Mazur.<cit.> The result can be written as

    η_RDBη_0,RDB = -2^1/3 + (3 Λ+√(2+9Λ^2))^2/32^2/3Λ(3Λ+√(2+9Λ ^2))^1/3.

[sZyLa]Λdimensionless shear rate (Sec. <ref> only)-first used in Eq. (<ref>)
[sZyLb]Λdimensionless elongation rate (Sec. <ref> only)-defined by Eq. (<ref>)
[sZyLc]Λdimensionless step-rate value (Sec. <ref> only)-defined by Eq. (<ref>)


 §.§ Discussion

The influence of salinity on the non-Newtonian viscosity of polyelectrolytes has been subject to extensive experimental investigations. Both a decrease in viscosity and a shift of the onset of shear-thinning towards higher shear rates with increasing salt concentration are well-known features of polyelectrolyte solutions. <cit.> It is also reported that polyelectrolytes containing more intrinsically rigid molecules demonstrate a larger shear-thinning. The C-FENE-P dumbbell model predicts these trends, as shown in Fig. <ref>(b).
The asymptotic value of the shear-thinning exponent (-2/3) is identical to that of the FENE-P dumbbell and FENE-P bead-spring-chain models<cit.> and matches experimental data for dilute partially hydrolyzed polyacrylamides. <cit.> One should note, however, that the slope of the log-log viscosity curve of the C-FENE-P dumbbells [see Fig. <ref>(b)] is changing gradually with the shear rate ranging from -2/3 to 0. The C-FENE-P model can therefore explain the shear-thinning exponent values in this range.
Much less data are available on the normal stress coefficients and their dependence on the salt concentration. The shape of the Ψ_1(γ̇) curve predicted by the C-FENE-P model is realistic and matches the qualitative description given by Bird et al.<cit.> The simple nonlinear relation between Ψ_1(γ̇) and η(γ̇) [Eq.  (<ref>)] has been tested experimentally for partially hydrolyzed polyacrylamides in a recent study by Lozhkina<cit.>; this relation was proven to be qualitatively correct for solutions of high-molecular-weight polyacrylamides but to not hold for their lower-molecular-weight counterparts.
The second normal stress coefficient vanishes in many kinetic theory-based polymer fluid models,<cit.> including the C-FENE-P dumbbells. In practice, it is reported that |Ψ_2(γ̇)|<<|Ψ_1(γ̇)|; hence, the second normal stress difference does not play a significant role for most flows of practical interest.<cit.>




§ STEADY EXTENSIONAL FLOW




 §.§ Analysis

The steady simple shearfree flow velocity field is

    v=[ -12ε̇ x_1 -12ε̇ x_2    ε̇ x_3 ],

[sZye1]ε̇elongation rates^-1first used in Eq. (<ref>)
where ε̇ is the time-independent elongation rate, which can take positive and negative values. Equation (<ref>) defines uniaxial extension at ε̇>0 and biaxial stretching at ε̇<0.
The rate-of-strain tensor, the stress tensor, and the Oldroyd derivative of the latter are all diagonal,

    γ̇    = diag (-ε̇,-ε̇,2ε̇), 
    τ    =diag (τ_11, τ_22, τ_33), 
    τ_(1)    = diag (τ_11, τ_22, -2τ_33 )ε̇,

with τ_11=τ_22 due to the flow symmetry. The only material function characterizing the fluid in this type of flow is the extensional viscosity. Following Bird et al.,<cit.> we denote it by η̅ and define by:

    τ_33-τ_11 = -η̅(ε̇)ε̇.

[sZyhg]η̅extensional viscosityPa·sdefined by Eq. (<ref>)
[sZyhh]η̅_0zero-elongation-rate extensional viscosityPa·sfirst used in Eq. (<ref>)
Substituting Eqs. (<ref>)-(<ref>) into the constitutive equation leads to

    b3Z τ_11 + λε̇τ_11   = nkT λε̇, 
    b3Z τ_33-2λε̇τ_33   = - 2nkTλε̇.

Having replaced τ_11, τ_33, and ε̇ with dimensionless quantities

    𝕋    = 2 τ_11+τ_33nkT, 
    𝔻   = τ_11-τ_33nkT, 
    Λ   = λε̇,
 
[sTa]𝕋dimensionless trace of the stress tensor-defined by Eq. (<ref>)
[xD]D_tmaterial derivatives^-1first used in Eq. (<ref>)
[xd]∂_tordinary time derivatives^-1first used in Eq. (<ref>)
[xZz]∇del operatorm^-1first used in Eq. (<ref>)
[sDb]𝔻dimensionless normal stress difference-defined by Eq. (<ref>)
one arrives after simple rearrangements at

    b3 Z 𝕋 + 2 Λ𝔻   = 0, 
    b3 Z 𝔻 + Λ (𝕋 - 𝔻)    = 3 Λ.

 
Equations (<ref>) and (<ref>) can be solved numerically. The impact of E on the extensional viscosity of the C-FENE-P dumbbells is illustrated in Fig. <ref>(b). It is seen that the C-FENE-P dumbbells show extensional thickening at ε̇>0. An increase in salinity (hence, a decrease in E) leads to an overall drop in the extensional viscosity. The local minimum in η̅ at negative elongation rates, which vanishes in the RDB limit, becomes more pronounced as salinity increases.
One can also keep track of the mean-square relative extension, x, which is shown in Fig. <ref>(a). The extension is smallest at equilibrium and follows the trends similar to those of η̅ at ε̇ >0.
Some properties of the extensional viscosity curves can be obtained by analytical means. In particular, Eq. (<ref>) can be used to calculate the zero-elongation-rate extensional viscosity, η̅_0. Considering the limit ε̇→ 0 and then applying Eqs. (<ref>) and (<ref>) yield

    η̅_0=3η_0,

as expected.
At very large positive or negative elongation rates (in the limit |ε̇|→∞), the relative extension of the dumbbells approaches one; hence, the asymptotic behavior of the C-FENE-P and FENE-P dumbbells must be identical to that of the RDB model. In the RDB limit, Eqs.  (<ref>) and (<ref>) become

    (1 - 𝕋3) 𝕋 + 2 Λ𝔻   = 0, 
     (1 - 𝕋3) 𝔻 + Λ (𝕋-𝔻)    = 3Λ.

This system has three solutions of which only one,

    𝕋    = 32(1-Λ-√(1-2Λ+9Λ^2)), 
    𝔻   = 34(-1+5Λ +√(1-2Λ+9Λ^2)),

provides 𝕋 = 𝔻 = 0 at equilibrium (Λ = 0) and hence is physically relevant. It is seen that 𝔻∼ 6 Λ at Λ→∞, and 𝔻∼ 3/2 Λ at Λ→ -∞. This means that for all the models under consideration (FENE-P, C-FENE-P, and RDB), η̅ approaches η̅_+∞=6nkTλ at large positive elongation rates and η̅_-∞=(3/2)nkTλ at large negative elongation rates. Equation (<ref>) also provides an exact analytical expression for the extensional viscosity in the RDB limit,

    η̅(Λ)η̅_0,RDB=-1+5Λ +√(1-2Λ+9Λ^2)4Λ.



 §.§ Discussion

The C-FENE-P model predicts monotonic extensional thickening at positive elongation rates. The S-shaped curves, obtained in earlier theoretical works of Dunlap and Leal<cit.> and Ait-Kadi et al.,<cit.> are not reproduced by our model. The impact of salinity on the extensional viscosity, as shown in Fig. <ref>, is in qualitative agreement with experimentally observed trends.<cit.> For the negative elongation rates (biaxial stretching), we found no experimental results for comparison.




§ SMALL-AMPLITUDE OSCILLATORY SHEAR FLOW




 §.§ Analysis

In small-amplitude oscillatory shear (SAOS) flow, the fluid  velocity field is given by

    v = 
    [ γ̇_ 12(t) x_2             0             0 ],

where the harmonically oscillating shear rate

    γ̇_ 12(t)=γ̇_ 21(t)=γ̇_ 0cos(ω t)

[sZyg01]γ̇_ 0shear rate amplitude (Sec. <ref> only)s^-1first used in Eq. (<ref>)
[sZyg02]γ̇_ 0step-rate value (Sec. <ref> only)s^-1first used in Eq. (<ref>)
[sZyz]ωangular frequencys^-1first used in Eq. (<ref>)
is the only independent non-zero component of the rate-of-strain tensor. Here, γ̇_ 0 is the amplitude of oscillations  and  ω is the angular frequency. The former is assumed very small so that the dependency of the stress tensor components on γ̇_ 0 is completely described by the lowest-order terms (first order for the shear stresses and second order for the normal stresses). The stress tensor of a fluid undergoing SAOS flow has the general form

    τ(t)=τ_11(t)    τ_12(t)     0 
    τ_12(t)    τ_22(t)     0 
    
    0     0    τ_33(t)
    ,

with 

    τ_(1)=∂_t τ-
    2τ_12(t)    τ_22(t)     0 
    τ_22(t)     0     0 
    
    0     0     0
    γ̇_ 12(t),

where ∂_t stands for ordinary time derivative.
For polymer solutions, the shear stress oscillates around zero with frequency ω; in contrast to Newtonian liquids, the oscillations are not in phase with those of the shear rate. The normal stress differences oscillate around generally non-zero mean values with the double frequency 2ω. <cit.> The properties of the fluid in SAOS flow are then described by eight material functions – η', η”; Ψ_1^d, Ψ'_1, Ψ”_1; and Ψ_2^d, Ψ'_2, Ψ”_2 – defined by

    τ_12   = -η'(ω) γ̇_ 0cos (ω t)-η”(ω) γ̇_ 0sin (ω t), 
    
    N_1    = - Ψ_1^d(ω) γ̇_ 0^ 2 - Ψ'_1(ω) γ̇_ 0^ 2cos (2ω t) - Ψ”_1(ω) γ̇_ 0^ 2sin(2ω t), 
    
    N_2    = - Ψ_2^d(ω) γ̇_ 0^ 2 - Ψ'_2 (ω)γ̇_ 0^ 2cos (2ω t) - Ψ”_2(ω) γ̇_ 0^ 2sin(2ω t).
 
[sZyhc]η'in-phase component of the complex viscosityPa·sdefined by Eq. (<ref>)
[sZyhd]η”out-of-phase component of the complex viscosityPa·sdefined by Eq. (<ref>)
[sZyV3]Ψ_1^dfirst normal stress displacement coefficientPa·s^2defined by Eq. (<ref>)
[sZyV4]Ψ'_1real part of the complex first normal stress coefficientPa·s^2defined by Eq. (<ref>)
[sZyV5]Ψ”_1imaginary part of the complex first normal stress coefficientPa·s^2defined by Eq. (<ref>)
[sZyV92]Ψ_2^dsecond normal stress displacement coefficientPa·s^2defined by Eq. (<ref>)
[sZyV93]Ψ'_2real part of the complex second normal stress coefficientPa·s^2defined by Eq. (<ref>)
[sZyV94]Ψ”_2imaginary part of the complex second normal stress coefficientPa·s^2defined by Eq. (<ref>)
The combined quantities η'+iη”, Ψ'_1+iΨ”_1, and Ψ'_2+iΨ”_2 are commonly referred to as complex viscosity, first normal stress coefficient, and second normal stress coefficient, respectively, while Ψ_1^d and Ψ_2^d are the normal stress displacement coefficients.
The in-phase ("real") component of complex viscosity, η', describes the direct response of the fluid and can be associated with energy loss due to dissipation. The out-of-phase ("imaginary") component, η”, arises because the long polymer molecules do not react instantly to rapid flow changes. This leads to a latency, described by a phase shift. This latency can be interpreted as "elasticity" of the flow and associated with energy storage. <cit.> The complex viscosity components, η' and η”, can also be replaced by the storage and loss moduli, G' and G”, defined, respectively, by

    G'    = η”ω, 
    
    G”   =η' ω.

[sGa]G'storage modulusPadefined by Eq. (<ref>)
[sGb]G”loss modulusPadefined by Eq. (<ref>)
These ("elastic") moduli are often measured in experiments.
The expressions for the SAOS flow material functions of the C-FENE-P dumbbells can be obtained analytically. Substituting Eqs. (<ref>) and (<ref>) into the constitutive equation (<ref>), keeping only the lowest-order terms in γ̇_ 12, combining and rearranging the scalar equations, one gets

    N_1 + λ_ e∂_t N_1-2 λ_ eτ_12γ̇_ 12=0, 
        N_2+ λ_ e∂_t N_2=0, 
       τ_12+λ_ e∂_t τ_12 = -η_0 γ̇_ 12,
 
[sZyl0e]λ_ e"experimental" time constantsdefined by Eq. (<ref>)
where we have introduced the "experimental" time constant λ_ e by

    λ_ e (b, E) = η_0nkT=λℱ(3/b,E/b),
 
with η_0 being the zero-shear-rate viscosity of the C-FENE-P dumbbells given by Eq. (<ref>). As shown in Fig. <ref>,  this time constant is a monotonically increasing function of E, ranging from 3λ/(b+3) at E=0 to λ at E→∞. 
Substituting Eqs. (<ref>)-(<ref>) into Eqs.  (<ref>)-(<ref>) results in identities, holding for all (ω t). This leads to eight algebraic equations for the SAOS flow material functions. Then, the latter can be written in scaled form as functions of the oscillarory-flow Deborah number,<cit.>

    De=λ_ eω.

The linear viscoelastic response is

    η'η_0   = 11+De^2, 
    η”η_0   = De1+De^2,  
    G'nkT   = De^21+De^2,
    G”nkT    = De1+De^2,

[sDe]DeSAOS Deborah number-defined by Eq. (<ref>)
which is identical to that of a Maxwell fluid with viscosity nkTλ_ e and time constant λ_ e. <cit.>
The material functions related to N_2 – Ψ_2^d, Ψ'_2, and Ψ”_2 – all vanish, while N_1 is described by

    Ψ_1^dΨ_1,0   = 12(1+De^2),
    Ψ'_1Ψ_1,0   = 1-2 De^22(1+De^2)(1+4De^2), 
    Ψ”_1Ψ_1,0   = 3De2(1+De^2)(1+4De^2),
 
[sZyV2]Ψ_1,0zero-shear-rate first normal stress coefficientPa·s^2first used in Eq. (<ref>)
where Ψ_1,0=2nkTλ_ e^2 is the zero-shear-rate first normal stress coefficient of the C-FENE-P dumbbells. 
The scaled material functions, given by Eqs. (<ref>)-(<ref>), are shown in Fig. <ref>, where they are plotted against De.



The first normal stress displacement coefficient, Ψ_1^d, is a decreasing function of frequency. At low frequencies, Ψ_1^d→Ψ_1,0/2, and at high frequencies, Ψ_1^d∼ω^-2.
The real component of the complex first normal stress coefficient, Ψ'_1, is the only SAOS material function of the C-FENE-P model, which can take negative values. At very low frequencies, Ψ'_1→Ψ_1,0/2. At low-to-moderate frequencies, Ψ'_1 decreases with frequency, becoming zero at De=1/√(2)≈0.707, and continues to decline until the minimum point,

    Ψ'_1 = (2√(2)3-1 )Ψ_1,0≈ -0.057Ψ_1,0,

is reached at 
    De = [3+√(2)2√(2)]^1/2≈ 1.249.

Thereafter, Ψ'_1 starts to increase with frequency, approaching zero from below, with Ψ'_1∼ω^-2 at high frequencies.
The imaginary component of the complex first normal stress coefficient, Ψ”_1, increases linearly at low frequencies (Ψ”_1∼ω), reaching a maximum value of 

    Ψ”_1 =
    9[3(√(73)-5)/2 ]^1/223+5√(73)Ψ_1,0≈ 0.316Ψ_1,0

at 

    De = 12[ √(73)-56]^1/2≈ 0.384,

and decays quickly at large frequencies (Ψ”_1∼ω^-3).
Expressions (<ref>)-(<ref>) show that the SAOS material functions of the C-FENE-P, FENE-P, and RDB fluid models can be written in the same form. However, this form is not suitable for visualizing the impact of E, since all the scaling factors depend on E. To investigate the E-dependence, we reformulate Eqs. (<ref>)-(<ref>) using salinity-independent scaling factors: the zero-shear-rate viscosity and the zero-shear-rate first normal stress coefficient of the FENE-P dumbbells, η_0,FENE and Ψ_1,0,FENE, in place of η_0 and Ψ_1,0, respectively, and λ in place of λ_ e.
The dependence of the in-phase complex viscosity component on E is shown in Fig. <ref>(a). It is seen that a decrease in E leads to a reduction in η' at lower frequencies but to an increase in η' at higher frequencies; and the onset of "frequency-thinning" is shifted towards higher frequency values as E decreases. 

The impact of E on the out-of-phase complex viscosity component is visualized in Fig. <ref>(b). A decrease in E mostly affects the low- and mid-frequency regions of the curves: the values of η” are reduced, and the maximum is shifted towards higher frequencies.  At very high frequencies, E has no effect on η”: the curves of the C-FENE-P, FENE-P, and RDB models are asymptotically identical.
The storage and loss moduli are found to depend on E in a simple way [see Figs. <ref>(c) and <ref>(d)]. If the curves are plotted using a log-log scale, a reduction in E results in a translation of both G' and G” curves to the right due to a decrease in the time constant; the maximum values of both moduli thus remain unchanged.

The E-dependence of the SAOS material functions related to N_1 is shown in Fig. <ref>. At low and moderate frequencies, Ψ_1^d is affected in a way similar to η [see Fig. <ref>(a)], while Ψ”_1 is affected in a way similar to η” [see Fig. <ref>(c)].
The impact of E on Ψ'_1 is more complex but still follows the same general trend: a decrease in E leads to an overall reduction in the magnitude of Ψ'_1 and shifts the characteristic points towards higher frequencies [see Fig. <ref>(b)]. At very high frequencies, all three material functions are insensitive to E.  


 §.§ Discussion

The qualitative shape and asymptotic behavior of the complex viscosity components η'(ω), η”(ω) and the elastic moduli G'(ω), G”(ω), as predicted by the C-FENE-P dumbbell model, are fully consistent with physical arguments provided by Bird et al.<cit.> A decrease in these material functions with increasing salinity at fixed ω, as shown in Fig. <ref>, is an experimentally observed feature of polyelectrolyte solutions.<cit.> 
Furthermore, the C-FENE-P model predicts a decrease of  the experimental time parameter, λ_ e, with increasing salinity (see Fig. <ref>); this leads to a shift of the characteristic points, occurring at fixed Deborah numbers, towards higher frequencies. This prediction is compatible with the experimental results of Ihebuzor,<cit.> who reported the G'-G” crossover frequency to increase with salinity.
At the same time, the C-FENE-P model is not capable of resolving the quantitative relations between the SAOS material functions, in particular, it predicts that the G'(ω) and G”(ω) curves intersect at the point of maximum of G”(ω), as follows from Eqs. (<ref>) and  (<ref>). This is not observed in experiments: as seen from the recent report of Ihebuzor,<cit.> G”(ω) continues to increase at frequency values larger than the crossover frequency. This mismatch, however, is to be expected: dumbbell models, with their single relaxation time, cannot properly describe the complex oscillatory motion of a real polymer molecule, with a wide spectrum of relaxation times, in cases when the flow pattern changes rapidly. An extension of the C-FENE-P dumbbell model to its bead-spring-chain variant might resolve this issue but would involve significant mathematical and computational complexity, which we are intentionally trying to avoid in this work.
Finally, we did not find any experimental works allowing for an analysis of the results related to the impact of salinity on the first normal stress difference in SAOS flow.




§ START-UP AND CESSATION OF STEADY SHEAR FLOW




 §.§ Analysis

Start-up and cessation of steady shear flow are two closely related transient shear flows. The velocity field in these flows is described by Eq.  (<ref>); the stress tensor and its Oldroyd derivative are given by Eqs. (<ref>) and (<ref>), respectively, while the only independent non-zero component of the rate-of-strain tensor is

    γ̇(t) ≡γ̇_ 12(t) = γ̇_ 21(t) = γ̇_ 0Θ(± t),
 
[xZi]ΘHeaviside step-function-first used in Eq. (<ref>)
where γ̇_ 0 is a constant and Θ(t) is the Heaviside step-function; for this reason, these flows are known in experimental rheology as step-rate tests. <cit.>
Choosing the positive sign in Eq. (<ref>) corresponds to the start-up case. The fluid is at rest at t<0, while a constant shear rate γ̇_ 0 is applied suddenly at t=0. After a while, shear and normal stresses build up and approach their steady-shear-flow values.
In contrast, the negative sign in Eq. (<ref>) yields the cessation case, which is the inverse situation. The fluid is flowing steadily with constant shear rate γ̇_ 0 at t<0 before the flow is instantaneously stopped (the shear rate is removed) at t=0. Shear and normal stresses decay as the fluid approaches equilibrium.
The material functions of the fluid in this kind of flow – η^±, Ψ^±_1, and Ψ^±_2 – are defined by

    τ_12    = -η^±(γ̇_ 0,t) γ̇_ 0, 
    
    N_1    = -Ψ^±_1(γ̇_ 0,t)γ̇_ 0^ 2,
    
    N_2    = -Ψ^±_2(γ̇_ 0,t)γ̇_ 0^ 2

[sZyhe]η^-shear stress relaxation functionPa·sdefined by Eq. (<ref>)
[sZyhf]η^+shear stress growth functionPa ·sdefined by Eq. (<ref>)
[sZyV6]Ψ_1^-first normal stress difference relaxation functionPa ·s^2defined by Eq. (<ref>)
[sZyV7]Ψ_1^+first normal stress difference growth functionPa ·s^2defined by Eq. (<ref>)
[sZyV95]Ψ_2^-second normal stress difference relaxation functionPa ·s^2defined by Eq. (<ref>)
[sZyV96]Ψ_2^+second normal stress difference growth functionPa ·s^2defined by Eq. (<ref>)
and are known as the shear stress, first normal stress, and second normal stress difference growth (+) or relaxation (-) functions, respectively. When presented graphically, these material functions are commonly normalized to their steady-state values η(γ̇_ 0), Ψ_1(γ̇_ 0), and Ψ_2(γ̇_ 0).
Substituting Eqs. (<ref>), (<ref>), and (<ref>) into the constitutive equation of the C-FENE-P dumbbells (<ref>) and introducing dimensionless quantities

    𝕋_ij   =τ_ijnkT, 
    
    r    = t/λ, 
    Λ   = λγ̇_ 0,

where i,j∈{1,2,3}, lead to a system of nonlinear ordinary differential equations:

    b3Z 𝕋_11+𝕋'_11-2 𝕋_12ΛΘ(± r) - (𝕋_11-1)Z'Z=0,  
       b3Z 𝕋_22+𝕋'_22 - (𝕋_22-1)Z'Z=0,   
       b3Z 𝕋_33+𝕋'_33 - (𝕋_33-1)Z'Z=0,   
       b3Z 𝕋_12 + 𝕋'_12 - 𝕋_22ΛΘ(± r)-𝕋_12Z'Z =  -ΛΘ(± r).

Here, primes denote differentiation with respect to r.

The expression for the Z-factor [Eq. (<ref>)] is rewritten accordingly to complete the system. 
The initial conditions are imposed at r=0, where the stress tensor components are set to zero (start-up case) or to their steady-shear-flow values (cessation case).
Since 𝕋_22=𝕋_33=0 both in equilibrium and in steady shear flow (see Sec. <ref>), it follows from Eqs. (<ref>) and (<ref>) that N_2(r)=0 identically. The rest of the system can be solved numerically. We have used Wolfram Mathematica for this purpose.
The results for the start-up case are shown in Fig. <ref>. At very low values of Λ, the material functions grow seemingly monotonically [see Figs. <ref>(a) and <ref>(b)]. At higher Λ, they typically undergo one or several oscillations around the steady-flow values before they stabilize; a stress overshoot, i.e.,  a time interval where the stresses are higher than their steady-flow values, is clearly seen in Figs. <ref>(c)-<ref>(f). Note that the overshoot is not only  observed for rigid dumbbells but also much more pronounced in the RDB limit. Hence, this phenomenon is not caused by dumbbell stretching, as one might have suggested. In addition, it is seen that at higher Λ, the relative magnitude of overshoots increases and the overshoots are shifted towards earlier times; it takes less time for the stresses to approach their steady-flow values, and for any fixed Λ, shear stress builds up and stabilizes faster than N_1.
The material functions describing the start-up case depend strongly on E. At low Λ, a decrease in E results in a general increase in η^+/η(γ̇_ 0) and Ψ^+_1/Ψ_1(γ̇_ 0) [see Figs. <ref>(a) and <ref>(b)].  At higher values of Λ, the impact of E on stress growth becomes more complex. For shear stress growth functions, an overshoot appears [see Figs <ref>(c) and <ref>(e)], and as E decreases, the overshoot is shifted towards later times (the higher the Λ, the more pronounced the effect). At the same time, the relative overshoot magnitude decreases with a decrease in E, reaching its minimal value in the FENE-P dumbbell limit (E=0). For N_1 growth, the situation is somewhat similar [see Figs. <ref>(d) and <ref>(f)]; however, the magnitude of N_1 overshoot can change non-monotonically with E: at high Λ, it increases, reaches a maximum, and then decreases as E decreases, as shown in Fig. <ref>(f).
In other words, an increase in solvent salinity can either increase or decrease η^+/η(γ̇_ 0), Ψ^+_1/Ψ_1(γ̇_ 0) and the relative N_1 overshoot, depending on values of Λ and E, but always decreases the relative shear stress overshoot.

The numerical results for the cessation case are presented in Fig. <ref>. Both the shear stress and the first normal stress difference relaxation functions decay monotonically, quickly approaching zero; the curves for η^-/η(γ̇_ 0) and Ψ_1^-/Ψ_1(γ̇_ 0) overlap. Equilibrium is approached faster at higher Λ, as seen from the comparison of Figs. <ref>(b) and <ref>(a).
At later times, the material functions decay exponentially. In Fig. <ref>, this is seen as regions where the curves become straight lines. The exponential decay is preceded by a region of faster decrease. An analysis of system (<ref>)-(<ref>) shows that the eigenvalue of its linearized version corresponding to 𝕋_12 equals -ℱ(3/b,E/b); hence, 

    η^-(t,γ̇_ 0)η(γ̇_ 0)∼Ψ_1^-(t,γ̇_ 0)Ψ_1(γ̇_ 0)∼exp[-ℱ(3/b,E/b)r] ∼exp (-t/λ_ e)

asymptotically at later times, where λ_ e is the time constant defined previously by Eq. (<ref>). Our numerical simulations confirm this analytical result.
The impact of E on η^- and Ψ_1^- is twofold. First, the region of fast decrease at early times, which is abrupt and step-like at large E, becomes smoother and less pronounced, as E decreases. Second, a decrease in E leads to a faster decay rate in the exponential regime at late times, since ℱ(3/b,E/b) is a decreasing function of E. Both effects are seen in Fig. <ref>.


 §.§ Discussion

At fixed salinity, the appearance of a shear stress overshoot at the start-up of steady shear flow and the way its magnitude and position depends on the step-rate value γ̇_ 0 are experimentally observed features of polymer solutions.<cit.>
Our results show that a decrease in salinity leads to a more pronounced shear stress overshoot, an effect that has been  observed experimentally by Zebrowski and Fuller.<cit.> 
In the case of steady shear flow cessation, the C-FENE-P model predicts a faster return to equilibrium at higher shear rates if salinity is held constant, as shown in Fig. <ref>. The same trend is observed in experiments.<cit.> Furthermore, our results show that at fixed shear rate value, an increase in salinity leads to a slower decay of shear stresses, which also maches the trend reported by  Zebrowski and Fuller.<cit.> 
Finally, Islam<cit.> investigated the shape of normalized shear stress relaxation functions for commercial partially hydrolyzed polyacrylamides. Two  clearly distinct regions – an exponential decay at later times preceded by a faster decrease at earlier times, as predicted by the C-FENE-P model – were observed. Furthermore, the decay rate in the exponential regime was found to be dictated by the polymer type and independent of γ̇_ 0 and polymer concentration, which is explained by our theoretical result expressed by Eq. (<ref>).
We are not aware of experimental data on the first normal stress difference growth and relaxation functions of polyelectrolyte solutions.




§ CONCLUSIONS


The predictions of the C-FENE-P dumbbell model for steady and transient shear and extensional flows are in a very good qualitative agreement with most of the experimentally observed trends for polyelectrolyte solutions. This involves the shape of the material functions, their dependence on flow parameters, how the material functions are affected by the solvent salinity, and the ways in which the material functions differ for more and less intrinsically rigid polyelectrolytes. At the same time, the C-FENE-P model is strikingly simple, being just slightly more mathematically complex than the original (uncharged) FENE-P dumbbell model. The closed-form constitutive equations allow us to investigate the rheology and fluid dynamics of steady and transient shear and extensional flows using a combination of analytical and simple numerical methods. In addition, the physical quantities of interest, such as the relative extension of the polyelectrolyte molecules, can be easily kept track of under analysis.
Concerning numerical simulation of complex flows, we assume that computer codes, which are capable of solving the equations of motion for a FENE-P dumbbell fluid, can be relatively easily adapted to tackle the constitutive equations of the C-FENE-P dumbbells. 
We believe that the C-FENE-P dumbbell model has all the potential to be a robust instrument, suitable for direct use in technological applications and helpful for qualitative understanding of flow phenomena in complex flows of dilute polyelectrolyte solutions.

This research has been funded by VISTA – a basic research program in collaboration between The Norwegian Academy of Science and Letters and Equinor. The authors acknowledge the Research Council of Norway and the industry partners (2019) of the National IOR Centre of Norway – ConocoPhillips Skandinavia AS, Aker BP ASA, Vår Energi AS, Equinor ASA, Neptune Energy Norge AS, Lundin Norway AS, Halliburton AS, Schlumberger Norge AS, Wintershall Norge AS, and DEA Norge AS – for support. D.S. is thankful to Tamara Shogina for numerous suggestions on improvement of this paper.




§ AVAILABILITY OF DATA

Data sharing is not applicable to this article as no new data were created or analyzed in this study.



