Viscoelastic Fluids

Viscoelastic fluids are materials that show both viscous behavior of a fluid and elastic behavior of a solid.

Viscoelastic fluids can exhibit qualitatively different phenomena from Newtonian fluid flow. These flow phenomena are due to the normal stresses that exist within viscoelastic fluids. Such well-known effects are:

Weissenberg effect or rod-climbing effect. Consider a rotating rod whose end is immersed into a viscoelastic fluid, for example bread dough. The rotation of the rod generates elastic forces and a consequent stretching of the fluid. This stretching results in a positive normal force and the fluid rises up the rod. By contrast, in Newtonian liquids inertial forces dominate and the fluid moves away from the rod. The Weissenberg effect is important in industrial mixing processes.
Die swell or extrudate swell. Consider extrusion of a polymeric liquid through a die. On exiting the die, the viscoelastic fluid swells, that is, it increases in shape and volume. This phenomenon is related to the normal stresses of the fluid inside the pipe. When the fluid exits the die, these normal stresses are relieved and the jet expands in the transverse direction.

Modeling viscoelastic fluids requires complex constitutive equations for the relation between the stress tensor and the velocity field. Simcenter STAR-CCM+ provides several non-linear viscoelastic models of the differential type. The constitutive equations for the tensor components form a system of non-linear partial differential equations that are discretized using the finite element method and solved using a direct solver, the viscous flow solver.

The viscous flow solver assumes fully developed velocity and stress at inlets and outlets. For an internal flow, fully developed velocity means that all velocity vectors are parallel to each other and perpendicular to the inlet and outlet cross sections. Fully developed stress means that the stress components are calculated based on the conditions considered for the fully developed velocity at inlet or outlet. With these conditions, the velocity and stress components in a pipe or channel with uniform cross section are independent of the streamwise direction. Notice that this last condition cannot be satisfied for three-dimensional flows of viscoelastic fluids with non-zero second normal stress difference due to the formation of secondary flow.

For 3D flows, the presence of the second normal stress brings about the secondary motion in the pipe. The streamlines are no longer perfectly parallel, the fluid gains velocity in the third direction, and there is no fully developed flow. As a result, fully developed conditions are not relevant for the three-dimensional flows of viscoelastic models with second normal stress. This remains true for three dimensional flows with the periodic condition as well.

The Oldroyd-B model, the Phan-Thien-Tanner model, and the Giesekus model can be used in various modes to fit the relaxation modulus accurately and to predict the shear thinning as well as normal stresses and extensional viscosity. The stresses of all of these models are split into a polymeric (viscoelastic) part $T^{p}$ and a Newtonian solvent part $T^{s}$ :

Figure 1. EQUATION_DISPLAY

T = T^{s} + T^{p}

(706)

where:

Figure 2. EQUATION_DISPLAY

T^{s} = 2 μ_{s} D, T^{p} = \sum_{i = 1}^{n} T_{i}

(707)

$μ_{s}$ is the solvent viscosity. $μ_{s}$ can be constant or dependent on the shear-rate. Any of the shear-rate dependent viscosity functions of Eqn. (701) to Eqn. (703) can be applied.
$D$ is the solvent viscosity rate-of-deformation tensor; see Eqn. (695).
$T_{i}$ is the stress tensor of the $i$ th mode as defined in the selected model.
$n$ is the number of modes.

To provide an estimate of the shear resistance of the viscoelastic material with respect to the local shear rate in the system, the viscoelastic equivalent viscosity is introduced:

Figure 3. EQUATION_DISPLAY

μ_{e q u i v} = \frac{{I I}_{T}}{{I I}_{D}} = \frac{{(t r [T])}^{2} - t r (T^{2})}{{(t r [D])}^{2} - t r (D^{2})}

(708)

The total viscosity is given by:

Figure 4. EQUATION_DISPLAY

μ_{t o t a l} {= μ}_{s} + μ_{e q u i v}

(709)

By default, the different models are expressed in terms of the stress tensor $T$ . For greater computational stability in simulations with high Weissenberg numbers, the models can be re-cast using the unique symmetric positive-definite square root $B$ of the conformation tensor $C = B^{T} \cdot B$ , related to $T$ by:

Figure 5. EQUATION_DISPLAY

T = G (B^{T} \cdot B - I)

(710)

where $G$ is the shear modulus and $I$ is the identity tensor. Theoretically, the conformation tensor $C$ is a symmetric positive-definite tensor. However, in the numerical simulation, the conformation tensor is not precisely positive definite, and in fact the onset of numerical instability at high values of Weissenberg numbers is linked to the loss of positive-definiteness of the conformation tensor [165]. The advantage of the square root formulation is the fact that the positive-definiteness of the conformation tensor is always guaranteed within the simulation and therefore it demonstrates substantial numerical stability at high Weissenberg number simulations.

This conformal form also uses the antisymmetric matrix $A$ that guarantees the square root conformation tensor $B$ remains symmetric positive-definite pointwise in space and time. The components of $A$ can be obtained according to the expressions given in the literature [150].

Oldroyd-B Model

The Oldroyd-B model [199] is mostly used to describe the rheological characteristics of polymer liquids composed at low concentration and moderate shear rates.

The Oldroyd-B model can be considered as a linear superposition of the Upper Convected Maxwell model stress for the polymeric part $T_{i}$ and a Newtonian solvent contribution $T^{s}$ .

Figure 6. EQUATION_DISPLAY

T_{i} + λ \overset{▿}{T_{i}} = 2 μ_{0} D

(711)

where:

$\overset{▿}{T}$ is the upper-convected derivative of the stress tensor $T$ :
Figure 7. EQUATION_DISPLAY

$\overset{▿}{T} = \frac{\partial}{\partial t} T + v \cdot \nabla T - {(\nabla v)}^{T} \cdot T - T \cdot \nabla v$

(712)
$μ_{0}$ is the zero-shear rate (model) dynamic viscosity.
$λ$ is the (modal) relaxation time.

The square-root conformal form of the Oldroyd-B model is:

Figure 8. EQUATION_DISPLAY

λ (\frac{\partial B}{\partial t} + v \cdot \nabla B - B \cdot \nabla v - A B) + \frac{1}{2} [B - {(B^{T})}^{- 1}] = 0

(713)

Linear Phan-Thien-Tanner Model

Figure 9. EQUATION_DISPLAY

[1 + \frac{ϵ λ}{μ_{0}} t r (T_{i})] T_{i} + λ [\overset{▿}{T_{i}} + ξ (T_{i} \cdot D + D \cdot T_{i})] = 2 μ_{0} D

(714)

where $ϵ$ is the Phan-Thien-Tanner parameter and $ξ$ is a model parameter that gives a combination of the upper-convected and the lower-convected derivatives of the stress tensor, $0 \leq ξ \leq 2$ .

When $ϵ$ is zero, the model reduces to the Johnson-Segalman model:

T_{i} + λ [\overset{▿}{T_{i}} + ξ (T_{i} \cdot D + D \cdot T_{i})] = 2 μ_{0} D

The square-root conformal form of the linear Phan-Thien-Tanner model is:

Figure 10. EQUATION_DISPLAY

λ (\frac{\partial B}{\partial t} + v \cdot \nabla B - B \cdot \nabla v - A B) + \frac{1 + ϵ (tr [B^{T} \cdot B] - 3)}{2} (B - {[B^{T}]}^{- 1}) = 0

(715)

The linear Phan-Thien-Tanner constitutive equation is compatible with the square-root conformal form only when $ξ = 0$ . If $ξ \neq 0$ , use the $T$ formulation.

Exponential Phan-Thien Tanner Model

Figure 11. EQUATION_DISPLAY

[\exp (\frac{ϵ λ}{μ_{0}} t r (T_{i}))] T_{i} + λ [\overset{▿}{T_{i}} + ξ (T_{i} \cdot D + D \cdot T_{i})] = 2 μ_{0} D

(716)

where $ϵ$ is the Phan-Thien-Tanner parameter and $ξ$ is a model parameter, as in the linear version of the model, but with an exponential factor in place of the first term. This model also reduces to the Johnson-Segalman model when $ϵ$ is zero.

The square-root conformal form of the exponential Phan-Thien-Tanner model is:

Figure 12. EQUATION_DISPLAY

λ (\frac{\partial B}{\partial t} + v \cdot \nabla B - B \cdot \nabla v - A B) + \frac{\exp [ϵ (tr [B^{T} \cdot B] - 3)]}{2} (B - {[B^{T}]}^{- 1}) = 0

(717)

The exponential Phan-Thien-Tanner constitutive equation is compatible with the square-root conformal form only when $ξ = 0$ . If $ξ \neq 0$ , use the $T$ formulation.

Giesekus-Leonov Model

Figure 13. EQUATION_DISPLAY

T_{i} + α \frac{λ}{μ_{0}} {T_{i}}^{2} + λ \overset{▿}{T_{i}} = 2 μ_{0} D

(718)

where $α$ is the Giesekus-Leonov parameter or mobility factor, $0 \leq α \leq 0.5$ .

The square-root conformal form of the Giesekus model is:

Figure 14. EQUATION_DISPLAY

λ (\frac{\partial B}{\partial t} + v \cdot \nabla B - B \cdot \nabla v - A B) + \frac{1}{2} (B - {[B^{T}]}^{- 1}) + \frac{a}{2} ((B - {[B^{T}]}^{- 1}) \cdot (B^{T} B - I)) = 0

(719)

eXtended Pom-Pom Model

The single-equation eXtended Pom-Pom model (XPP) is used for flows of branched polymers in complex cases involving both extension and shear flow. It cannot be derived from the general model.

Figure 15. EQUATION_DISPLAY

f (Λ, T_{i}) T_{i} + λ_{b} \overset{▿}{T_{i}} + G (f (Λ, T_{i}) - 1) I + \frac{α}{G} (T_{i} \cdot T_{i}) = 2 μ_{0} D

(720)

The function $f (Λ, T_{i})$ is given by

Figure 16. EQUATION_DISPLAY

f (Λ, T_{i}) = \frac{2}{ϵ} e^{\frac{2}{Q} (Λ - 1)} (1 - \frac{1}{Λ^{2}}) + \frac{1}{Λ^{2}} (1 - \frac{α tr (T_{i} \cdot T_{i})}{3 G^{2}})

(721)

where:

$Λ = \sqrt{1 + \frac{t r (T_{i})}{3 G}}$ is the stretch factor for the material.
$ϵ = \frac{λ_{s}}{λ_{b}}$ , where $λ_{s}$ and $λ_{b}$ are the relaxation times for the backbone stretch and backbone tube orientation of the polymer.
$G$ is the shear modulus.
$α$ is the anisotropy parameter.
$Q$ is the number of arms on the polymer molecule.

The square-root conformal form of the XPP model is:

Figure 17. EQUATION_DISPLAY

λ (\frac{\partial B}{\partial t} + v \cdot \nabla B - B \cdot \nabla v - A B) + \frac{1}{2} (f (C) B - {[B^{T}]}^{- 1}) + \frac{a}{2} ((B - {[B^{T}]}^{- 1}) \cdot (B^{T} B - I)) = 0

(722)

where:

f (C) = \frac{2}{ϵ} e^{\frac{2}{Q} (Λ - 1)} (1 - \frac{1}{Λ^{2}}) + \frac{1}{Λ^{2}} (1 - \frac{α t r {(B^{T} B - I)}^{2}}{3})

and

Λ = \sqrt{\frac{t r (B^{T} B)}{3}}

Rolie-Poly Model

The Rolie-Poly (ROuse LInear Entangled POLYmers) model, developed by Likhtman and Graham [188], is used for inhomogeneous flows of entangled polymers. The model is suitable for prediction of linear polymers such as polystyrene or LLDPE in complex flow fields including both shear and elongational flows.

The Rolie-Poly equation is:

Figure 18. EQUATION_DISPLAY

λ \overset{▿}{T_{i}} + [1 + \frac{λ f_{r e t r}}{λ_{R}} (1 + f_{c c r})] T_{i} + \frac{μ_{0} f_{r e t r}}{λ_{R}} I = 2 μ_{0} D

(723)

where:

$λ$ is the relaxation time.
$λ_{R}$ is the Rouse relaxation time.
$f_{r e t r}$ and $f_{c c r}$ are terms accounting for chain retraction and convective constraint, respectively:

$f_{r e t r}$ and $f_{c c r}$ are given by:

Figure 19. EQUATION_DISPLAY

\begin{matrix} f_{r e t r} = 2 (1 - \frac{1}{\sqrt{\frac{tr (T_{i})}{3 G} + 1}}) \\ f_{c c r} = β {(\frac{tr (T_{i})}{3 G} + 1)}^{δ} \end{matrix}

(724)

where $δ$ and $β$ are parameters measuring the convective constraint release mechanism. Note that $δ$ is negative and often considered a constant, -0.5, and that $β \geq 0$ .

Likhman and Graham [188] showed that, in the limit as $λ_{R} \to 0$ , the Rolie-Poly constitutive equation can be simplified to a non-stretching limit:

Figure 20. EQUATION_DISPLAY

λ \overset{▿}{T_{i}} + T_{i} + \frac{2 λ}{3} [t r (\nabla v) + \frac{1}{G} t r (\nabla v \cdot T_{i})] [G I + (1 + β) T_{i}] = 2 μ_{0} D

(725)

When $λ_{R} = 0$ , this equation is used instead of Eqn. (723).

The square-root conformal form of the Rolie-Poly model is:

Figure 21. EQUATION_DISPLAY

λ (\frac{\partial B}{\partial t} + v \cdot \nabla B - B \cdot \nabla v - A B) + \frac{1}{2} (B - {[B^{T}]}^{- 1}) + \frac{f_{retr} λ}{2 λ_{R}} (B + f_{ccr} (B - {[B^{T}]}^{- 1})) = 0

(726)

where

f_{retr} = 2 (1 - \sqrt{\frac{3}{tr (B^{T} \cdot B)}})

and

f_{ccr} = β {(\frac{tr (B^{T} \cdot B)}{3})}^{δ}

The square-root conformal form of the non-stretching Rolie-Poly model is:

λ (\frac{\partial B}{\partial t} + v \cdot \nabla B - B \cdot \nabla v - A B) + \frac{1}{2} (B - {[B^{T}]}^{- 1}) + \frac{λ tr [\nabla v \cdot (B^{T} \cdot B)]}{3} (B + β (B - {[B^{T}]}^{- 1})) = 0

(727)