Considering a Non-Constant Anisotropicity Parameter in the Giesekus Model

1. Introduction

The use of plastic in almost every aspect of our everyday lives positions the plastic industry as one of the most important industries. Thus, the optimization of processes producing polymeric materials is of the utmost importance to avoid process instabilities [1,2]. Such a task necessitates the execution of numerical simulations that at heart employ an accurate and reliable constitutive rheological model that is based on molecular principles [2].The Giesekus model is widely regarded as one of the most successful constitutive equations for polymer solutions and melts, as it can predict a nonzero second normal stress difference, stress overshoot at the onset of simple shear flow, and finite extensional viscosity [3,4,5,6,7]. It introduces the notion that in sufficiently concentrated solutions and polymer melts, the frictional properties should be anisotropic [8,9], controlled by the magnitude of the anisotropicity or Giesekus parameter, α. By using the phase-space theory for concentrated polymer solutions and polymer melts of Curtiss and Bird [10,11,12], Giesekus [13] showed that his anisotropicity parameter can be related to the so-called link tension coefficient ε via $\alpha =\left(1-\epsilon \right)/\left(2+\epsilon \right)$. Stephanou et al. [14,15] provided the theoretical basis indicating that the link tension coefficient should not be considered a constant but should depend on the nematic order parameter. The revised framework predicts a vanishing ε at equilibrium, which remedies the erroneous predictions of the Curtiss-Bird model when considering a constant ε, namely the approach of the transient shear viscosity as t→0 to a constant value and the spurious sign changes of the transient second normal stress coefficient in startup simple shear flow.

Given the modification proposed by Stephanou et al. [14,15], in this work, we herein revisit the original Giesekus model by considering a non-constant anisotropicity parameter following the relationship between α and ε, and the relation between the former and the nematic order parameter as proposed by Stephanou et al. [14,15]. As such, at equilibrium, the anisotropicity parameter is equal to ½. We elaborate on the implications on the modified Giesekus predictions when considering a variable anisotropicity parameter and discuss whether its predictions are more in line with rheological experimental data. In Sec. 2, we provide the original and modified Giesekus constitutive model, and proceed to Sec. 3, where we present a comparison between their predictions. We conclude in Sec. 4 where the most important outcomes of our work are summarized and where we provide our future plans.

2. Constitutive Model

We consider an isothermal and incompressible flow. The Giesekus model in this case reads as follows [16]:

$${\dot{c}}_{\alpha \beta ,[1]}=-{\displaystyle \frac{1}{{\tau }_R}}\left[{\displaystyle \frac{1}{2}}\left(c_{\alpha \gamma }{\beta }_{\beta \gamma }+c_{\beta \gamma }{\beta }_{\alpha \gamma }\right)-{\beta }_{\alpha \beta }\right].$$

(1a)

where c is the dimensionless conformation tensor and we have defined the upper-convected time derivative:

$${\dot{c}}_{\alpha \beta ,[1]}\equiv {\displaystyle \frac{\partial c_{\alpha \beta }}{\partial t}}+u_{\gamma }{\nabla }_{\gamma }{c}_{\alpha \beta }-c_{\alpha \gamma }{\nabla }_{\gamma }{u}_{\beta }-c_{\gamma \beta }{\nabla }_{\gamma }{u}_{\alpha }.$$

(1b)

For the dimensionless mobility tensor, we consider the linear expression considered by Giesekus [13,16], namely,

$${\beta }_{\alpha \beta }={\delta }_{\alpha \beta }+\alpha \left(c_{\alpha \beta }-{\delta }_{\alpha \beta }\right).$$

(1c)

where $\alpha$ is the anisotropicity or Giesekus parameter. Also, the stress tensor is given as,

$${\sigma }_{\alpha \beta }=G\left(c_{\alpha \beta }-{\delta }_{\alpha \beta }\right).$$

(2)

In this work, which we believe to be the first such work,we consider the Giesekus parameter to be given explicitly by the expression relating it to the link tension coefficient,

$$\begin{array}{l}\alpha ={\displaystyle \frac{1-\epsilon }{2+\epsilon }}\\ \epsilon ={\epsilon }_0{S}_2^2\end{array},$$

(3a)

where ${\epsilon }_0$ is a constant coefficient and

$$S_2^2={\displaystyle \frac{3}{2}}\left[\mathrm{tr}\left({〈uu〉}_{\mathrm{an}}\cdot {〈uu〉}_{\mathrm{an}}\right)\right],$$

(3b)

is the nematic order parameter [17] which depends on the anisotropic orientation tensor ${〈uu〉}_{\mathrm{an}}$ defined as [14,15]

$${〈uu〉}_{\mathrm{an}}\equiv 〈uu〉-{\displaystyle \frac{1}{3}}I.$$

(3c)

where I is the unit tensor and u the unit end-to-end vector of polymer chains. Note that the link tension coefficient is not constant: it vanishes at equilibrium, it increases with time in startup flows, and its steady-state value increases with the imposed flow rate [14]. By considering $〈uu〉=c/\mathrm{tr}c$, then

$$S_2^2={\displaystyle \frac{3}{2}}\left[{\displaystyle \frac{\mathrm{tr}{c}^2}{{\left(\mathrm{tr}c\right)}^2}}-{\displaystyle \frac{1}{3}}\right].$$

(3d)

As a final note, since the order parameter asymptotes at large shear rate to unity, the Giesekus parameter approaches the value $\alpha =\left(1-{\epsilon }_0\right)/\left(2+{\epsilon }_0\right)$; since thermodynamic admissibility dictates that the Giesekus parameter should be $0\le \alpha \le 1$ [16,18,19,20] we need to bound ${\epsilon }_0\le 1$. Finally, note that when α=0, the Giesekus model simplifies to the Upper-Convected Maxwell (UCM) one.

Asymptotic Behavior of the Model for Steady-State and Transient Shear Flow

We then analyze the asymptotic behavior of the two versions of the Giesekus model in the case of a steady and startup simple shear flow (SSF), described by the kinematics $u=\left(\dot{\gamma }y,0,0\right)$ where y denotes the Cartesian coordinate and $\dot{\gamma }$ the shear rate. We present the following material functions: the shear viscosity $\eta ={\sigma }_{\mathrm{yx}}/\dot{\gamma }$ and the first and second normal stress coefficients,${\Psi }_1=\left({\sigma }_{\mathrm{xx}}-{\sigma }_{\mathrm{yy}}\right)/{\dot{\gamma }}^2$ and ${\Psi }_2=\left({\sigma }_{\mathrm{yy}}-{\sigma }_{\mathrm{zz}}\right)/{\dot{\gamma }}^2$, respectively. The asymptotic behaviour is the same as the original Giesekus model, minus the fact that in the modified model that at equilibrium α=1/2 irrespective of the choice of ${\epsilon }_0$. Thus, for the original model

$$\begin{array}{l}{c}_{xx}=1+\left(2-\alpha \right){\mathrm{Wi}}^2\\ c_{xy}=\mathrm{Wi}\\ c_{yy}=1-\alpha {\mathrm{Wi}}^2\\ c_{zz}=1\end{array},$$

(4a)

whereas for the modified model,

$$\begin{array}{l}{c}_{xx}=1+{\displaystyle \frac{3}{2}}{\mathrm{Wi}}^2\\ c_{xy}=\mathrm{Wi}\\ c_{yy}=1-{\displaystyle \frac{1}{2}}{\mathrm{Wi}}^2\\ c_{zz}=1\end{array}.$$

(4b)

We thus obtain the zero-shear-rate material functions: for the original model

$$\begin{array}{l}{\eta }_0=G{\tau }_R\\ {\Psi }_{1,0}=2{\eta }_0{\tau }_R\\ -{\Psi }_{2,0}=\alpha {\eta }_0{\tau }_R\end{array},$$

(5a)

and for the modified model,

$$\begin{array}{l}{\eta }_0=G{\tau }_R\\ {\Psi }_{1,0}=2{\eta }_0{\tau }_R\\ -{\Psi }_{2,0}={\displaystyle \frac{1}{2}}{\eta }_0{\tau }_R\end{array}.$$

(5b)

Upon inception of the simple shear flow, we follow the methodology of Stephanou et al. [20] to obtain the explicit solutions for the time-dependent viscometric functions in the linear viscoelastic (LVE) limit. The results are the same as those of the original Giesekus model,

$$\begin{array}{l}{\eta }^{+}(t)={\eta }_0\left[1-\exp \left(-{\displaystyle \frac{t}{{\tau }_R}}\right)\right]\\ {\Psi }_1^{+}(t)={\Psi }_{1,0}\left[1-\left(1+{\displaystyle \frac{t}{{\tau }_R}}\right)\exp \left(-{\displaystyle \frac{t}{{\tau }_R}}\right)\right]\\ {\Psi }_2^{+}(t)=\alpha {\eta }_0{\tau }_R\left[1-{\displaystyle \frac{2t}{{\tau }_R}}\exp \left(-{\displaystyle \frac{t}{{\tau }_R}}\right)-\exp \left(-{\displaystyle \frac{2t}{{\tau }_R}}\right)\right]\end{array},$$

(6a)

except for the fact that the LVE envelope of the second normal stress coefficient for the modified version differs,

$${\Psi }_2^{+}(t)={\displaystyle \frac{1}{2}}{\eta }_0{\tau }_R\left[1-{\displaystyle \frac{2t}{{\tau }_R}}\exp \left(-{\displaystyle \frac{t}{{\tau }_R}}\right)-\exp \left(-{\displaystyle \frac{2t}{{\tau }_R}}\right)\right].$$

(6b)

3. Results and Discussion

3.1. Model Predictions in Steady-State Shear Flow

Figure 1 illustrates the variations in the conformation tensor components as functions of the dimensionless shear rate (Wi) in the case of steady-state SSF when the original Giesekus model is compared with the modified model at two different values of $\epsilon$ or ${\epsilon }_0$ (0.1 and 0.7). The behavior of $c_{xx}$as a function of Wi is the same for both models when small values of $\epsilon$ or ${\epsilon }_0$ (=0.1) are used. Similarly, as expected for small Wi (Wi < 0.1), $c_{xx}$ remains close to its equilibrium value of unity in both models and then increases as Wi increases. This is expected since in the linear regime the mobility tensor is equal to the unit tensor. However, when $\epsilon$ or ${\epsilon }_0$ is raised to 0.7, we note that the two models provide different predictions in the intermediate range (0.1 < Wi < 100) with the modified model consistently predicting lower $c_{xx}$values than the original model. However, at large shear rates, the two models agree again, since in this case $S_2^2\approx 1$ and due to Eq. (3a) $\epsilon \approx {\epsilon }_0$. A similar behavior is noted also for $c_{xy}$ [Figure 1(b)] and $c_{yy}$[Figure 1(c)]: in both cases, the predictions from the two models are the same at small and large shear rates, whereas the prediction of the modified model is below the original model’s one in the intermediate range. The scaling at large shear rates is the same for both models: $c_{xx}$ increases as ${\mathrm{Wi}}^{\frac{1}{2}}$ , $c_{xx}$ reaches an asymptotic value, and $c_{yy}$ decreases as ${\mathrm{Wi}}^{-\frac{1}{2}}$. Then, Figure 2 shows the variations of the material viscometric functions as a function of Wi for both the modified and original models. For the shear viscosity and the first normal stress coefficient, we note a similar behavior as for the conformation tensor: the two models provide similar predictions for $\eta$ and ${\Psi }_1$ for small and large shear rates, but the prediction of the modified model is below that of the original model predictions in the intermediate range. On the contrary, the behavior in the scaled second normal stress coefficient, $-{\Psi }_2/\left({\eta }_0{\tau }_R\right)$, is quite different [Figure 2(c)]: the two models provide a similar prediction for a small value of $\epsilon$ or ${\epsilon }_0$ (= 0.1), with a similar zero-shear rate value (equal to 9/21 for the original model and ½ for the modified model) as dictated by Eqs. (5), and the scaling at large shear rates is $\eta \propto {\mathrm{Wi}}^{-1}$, ${\Psi }_1\propto {\mathrm{Wi}}^{-\frac{3}{2}}$, and ${\Psi }_2\propto {\mathrm{Wi}}^{-2}$ (Table 7.3-5 of Bird et al. [6]); however, they differ significantly at low and intermediate shear rates at larger values of $\epsilon$ or ${\epsilon }_0$ (= 0.7): the original Giesekus model predicts $-{\Psi }_2/\left({\eta }_0{\tau }_R\right)$ to be equal to 1/9, whereas the modified version always predicts it to be equal to ½. This also has implications to the startup behavior as will be described below.

3.2. Model Predictions in Steady-State Uniaxial Elongational Flow

In the case of uniaxial elongational flow (UEF), described by the kinematics $u=\left(\dot{\epsilon }x,-\frac{1}{2}\dot{\epsilon }y,-\frac{1}{2}\dot{\epsilon }z\right)$, where $\dot{\epsilon }$ is the elongation rate (x is the flow direction, y is the velocity gradient direction, and z is the neutral direction), the material function to be analyzed is the elongational viscosity, ${\eta }_{\mathrm{E}}=\left({\sigma }_{\mathrm{xx}}-{\sigma }_{\mathrm{yy}}\right)/\dot{\epsilon }$. A simple analytical solution can easily be obtained:

$$\begin{array}{l}{c}_{xx}={\displaystyle \frac{1}{2\alpha }}\left\{\sqrt{1+4{\mathrm{Wi}}^2-4\alpha \left(1-2\alpha \right)\mathrm{Wi}}-\left(1-2\alpha -2\mathrm{Wi}\right)\right\}\\ c_{yy}={\displaystyle \frac{1}{2\alpha }}\left\{\sqrt{1+{\mathrm{Wi}}^2+4\alpha \left(1-2\alpha \right)\mathrm{Wi}}-\left(1-2\alpha +\mathrm{Wi}\right)\right\}\end{array}.$$

(7a)

so that

$${\displaystyle \frac{{\eta }_{\mathrm{E}}}{{\eta }_0}}={\displaystyle \frac{1}{2\alpha }}\left\{3+{\displaystyle \frac{1}{\mathrm{Wi}}}\sqrt{1+4{\mathrm{Wi}}^2-4\alpha \left(1-2\alpha \right)\mathrm{Wi}}-\sqrt{1+{\mathrm{Wi}}^2+4\alpha \left(1-2\alpha \right)\mathrm{Wi}}\right\}.$$

(7b)

The result Eq. (7b) is well-known (see Table 7.3-5 with λ₂=0 in Bird et al. [6]). In the case of the original Giesekus model, the Giesekus parameter is a constant, whereas in the modified version it is given from Eq. (3a). In both models, the elongation viscosity reaches the asymptote ${\eta }_E/{\eta }_0=2/\alpha$ as we may easily note from Eq. (7b) (also see Eq. (7B.8-4) of Bird et al. [6]).

Figure 3 shows the conformation tensor elements $c_{xx}$, $c_{yy}$, and the dimensionless elongation viscosity (${\eta }_E/{\eta }_0$) versus Wi. In Figure 3(a), $c_{xx}$ starts from its equilibrium value of unity at low Wi for both models, and then increases with Wi, with higher values as $\epsilon$ or ${\epsilon }_0$ increases to 0.7. Similarly, in Figure 3(b), $c_{yy}$ also starts from its equilibrium value of unity at low Wi for both models, but then decreases as the shear rate increases, exhibiting higher values as $\epsilon$ or ${\epsilon }_0$ increases. The elongation viscosity starts from the value of 3, as per Trouton’s law, and increases as the elongation rate increases, exhibiting smaller values as $\epsilon$ or ${\epsilon }_0$ increases [Figure 3(c)]. When ${\epsilon }_0$= 0.1, both models provide similar predictions for both $c_{xx}$and $c_{yy}$, whereas the prediction of the modified model is below the prediction of the original model for $c_{xx}$; on the contrary, it seems that the prediction for $c_{yy}$ is the same for both models, even with a larger value of $\epsilon$ or ${\epsilon }_0$. The predictions for ${\eta }_E/{\eta }_0$ are noted to be the same for both models at large values of $\epsilon$ or ${\epsilon }_0$ whereas at smaller values, the modified version is below the original Giesekus model’s prediction in the intermediate range of elongation rates. The scaling at large shear rates is the same for both models, namely ${\tilde{c}}_{xx}\propto \mathrm{Wi}$ and ${\tilde{c}}_{yy}\propto {\mathrm{Wi}}^{-1}$, whereas, as mentioned above, the elongation viscosity always reaches the asymptote ${\eta }_E/{\eta }_0=2/\alpha$.

3.3. Model Predictions in Start-Up Shear Flow

Figure 4 and Figure 5 show the dimensionless material functions ${\eta }^{+}/{\eta }_0$, ${\Psi }_1^{+}/\left({\eta }_0{\tau }_R\right)$, and $-{\Psi }_2^{+}/\left({\eta }_0{\tau }_R\right)$ for startup SSF as a function of dimensionless time for ${\epsilon }_0$ = 0.1 and ${\epsilon }_0$ = 0.7, respectively, at Wi = 1, 10, and 100. In Figures 4(a) and 4(b), the predictions of ${\eta }^{+}/{\eta }_0$ and ${\Psi }_1^{+}/\left({\eta }_0{\tau }_R\right)$ for both models exhibit similar overall trends for all Wi: they increase at early times following the LVE envelope, Eqs. (6), reach an overshoot and then decrease, before reaching their steady state value. Also, both models' predictions the shear viscosity present an undershoot after the overshoot, which is more easily noticeable at Wi=100. This oscillatory behavior for the original Giesekus model was first noted by Giesekus [8]. Such undershoots, first provided by the experimental work of Costanzo et al. [21], has been associated with the tumbling of polymer chains in simple shear [22,23]. When $\epsilon$ or ${\epsilon }_0$ is small (=0.1), the predictions of the original and modified Giesekus models are identical, as seen in Figure 4, whereas at a larger value (=0.7), there are clear differences between the original and modified models (Figure 5). In particular, the predictions of the modified Giesekus model for ${\eta }^{+}/{\eta }_0$ and ${\Psi }_1^{+}/\left({\eta }_0{\tau }_R\right)$ are noted to be below the original Giesekus model predictions at large times for intermediate shear rates (Wi=1 and 10), whereas the predictions are similar at larger shear rates irrespective of time. On the other hand, the predictions of the two models for $-{\Psi }_2^{+}/\left({\eta }_0{\tau }_R\right)$ are noted to be completely different even in the LVE regime, as indicated by the different expressions for the LVE envelope [see Eqs. (6)]. The modified model is seen to shift upwards following Eq. (6b) and is noted to predict a larger value than the original model at small shear rates (Wi=1). However, at larger shear rates (Wi=10 and 100), both models exhibit the same steady-state value. This behavior is more in line with experimental data [14]; as seen in Figure 8(b) of Stephanou et al. [20], the LVE envelope predicted by the original Giesekus model (note that the LVE envelope provided there is not exactly the same as the original Giesekus model provided here; however, due to the very small value of the slip parameter considered in their work, ξ=0.03, their LVE envelope should not differ substantially from the original Giesekus model one) is below the LVE envelope suggested by the experimental data. Thus, the use of a variable Giesekus parameter, as done in the modified model, would certainly improve the comparison against said experimental data.

4. Conclusions

The Giesekus model has accumulated since its introduction in the early 1980s increased attention due to its simplicity and its capacity to conform to various experimental attributes [3]. Giesekus [13], to validate his postulate, showed that his expression for the mobility tensor can be obtained through a linearization of the mobility tensor proposed in phase-space theory for concentrated polymer solutions and polymer melts of Curtiss and Bird [10,11,12], which allowed him to relate his anisotropicity parameter to the link tension coefficient of Curtiss-Bird. However, recent evidence by Stephanou et al. [14,15] unambiguously shows that the link tension coefficient should not be considered a constant but should depend on strain rate and time via the nematic order parameter. Since the implications of using a variable anisotropicity parameter have not been presented before, in this work, we provide the predictions of a modified Giesekus parameter where the Giesekus parameter is no longer constant. Under steady-state simple shear flow, we note that the two models exhibit the same behaviour when $\epsilon$ or ${\epsilon }_0$ is small, and at any value of $\epsilon$ or ${\epsilon }_0$ at small and large shear rates. A similar behavior is noted for the transient shear viscosity and the transient first normal stress coefficient in start-up simple shear flow; however the predictions are noted to be very different in the case of the transient second normal stress coefficient in start-up simple shear flow where the LVE envelope is noted to shift to larger values, which is more in line with experimental data [14,20]. We thus expect that the use of a variable Giesekus parameter will improve the predictive capacity of constitutive models that employ the Giesekus parameter to more accurately predict the rheological behaviour of polymeric systems. We aim to undertake such an improvement of the Stephanou et al. [20] in the future.