Analysis of the Shear-Thinning Viscosity Behavior of the Johnson–Segalman Viscoelastic Fluids

Abstract

This paper presents a numerical comparison of viscoelastic shear-thinning fluid flow using a generalized Oldroyd-B model and Johnson–Segalman model under various settings. Results for the standard shear-thinning generalization of Oldroyd-B model are used as a reference for comparison with those obtained for the same flow cases using Johnson–Segalman model that has specific adjustment of convected derivative to assure shear-thinning behavior. The modeling strategy is first briefly described, pointing out the main differences between the generalized Oldroyd-B model (using the Cross model for shear-thinning viscosity) and the Johnson–Segalman model operating in shear-thinning regime. Then, both models are used for blood flow simulation in an idealized stenosed axisymmetric vessel under different flow rates for various model parameters. The simulations are performed using an in-house numerical code based on finite-volume discretization. The obtained results are mutually compared and discussed in detail, focusing on the qualitative assessment of the most distinct flow field differences. It is shown that despite all models sharing the same asymptotic viscosities, the behavior of the Johnson–Segalman model can be (depending on flow regime) quite different from the predictions of the generalized Oldroyd-B model.

1. Introduction

Many fluids of practical interest exhibit a complex behavior that cannot be predicted using mathematical models employing the classical Newtonian rheological laws. Phenomena such as shear-thinning/thickening, yield stress, stress relaxation or viscoelastic behavior are quite commonly observed in real fluids, but fail to be properly represented using classical Newtonian fluids models. A wide class of so-called non-Newtonain models was developed and used to capture specific fluid properties and flow behavior. A comprehensive overview and discussion of complex fluid rheology and corresponding models can be found for example in classical books [1,2,3] or in papers [4,5]. From the plethora of non-Newtonian fluids properties we will only remind and discuss two, the shear-thinning and viscoelasticity, that are relevant within the scope of this paper.

The shear-thinning behavior is typically captured by a specific sub-class of the so called generalized Newtonian models. In classical Newtonian models the stress tensor is directly proportional (by a constant coefficient named viscosity) to the fluid rate of strain tensor (which is nothing but symmetric part of velocity gradient). The generalized Newtonian models follow this concept, but allow the proportionality coefficient (viscosity) to be variable, typically depending on some relevant physical quantities, most importantly the invariants of the rate of strain tensor. Classical representative of this class of generalized Newtonian models is the well known power law viscosity model. The shear-thinning models then form a special sub-class where the apparent shear viscosity decays with increasing shear rate. Such kind of non-linear viscosity behavior is typical for many fluids in biomedicine (blood for example), food processing (ketchup, cream) or industry (some paints). The variable apparent viscosity is thus used and understood as a bulk parameter characterizing the complex microstructure and local response of the fluid.

Some effects in fluids, such as extrudate swell, rod climbing or Weissenberg effect, that are observed experimentally (see the monograph [6]), cannot be explained or predicted using Newtonian models (not even generalized Newtonian models). The reason is that in these situations the stress tensor is no more directly proportional to the rate of strain tensor. Instead, the fluid exhibits a complex viscoelastic behavior, merging the viscous and elastic response to applied strain history. Many suitable models have been developed to properly capture and predict the behavior of viscoelastic fluids. Various rate type or integral models are described, e.g., in [4,7,8]. The simplest from viscoelastic fluid rate type models are the Maxwell and Oldroyd type models [9]. They have fixed (constant) coefficients and besides of the standard viscous effects they can also describe the stress relaxation. Such viscoelastic behavior is again quite common in food industry, polymer processing (polymer injection molding) or biomedicine (synovial fluids or blood).

In many fluids several non-Newtonian properties (such as shear-thinning and viscoelasticity) are observed simultaneously, being more or less prominent at different flow regimes. This is typical for example in the so called polymeric fluids or in blood. Whole blood exhibits complex rheological behavior, namely shear-thinning viscosity ([10,11,12]) and deformation dependent viscoelastic behavior, at low shear rates, e.g., [13], or normal stress differences [14].

A number of nonlinear constitutive equations of differential and rate type have been considered in the literature to model blood in the vascular system as a shear-thinning and viscoelastic fluid, but none of them captures its complex behavior in a single framework. For example, the empirical five-constant generalization of the Oldroyd-B model introduced by Yeleswarapu et al. in [15] has been obtained by fitting experimental data to one-dimensional flow. In that model the constant viscosity in the Oldroyd-B model is replaced by a generalized Newtonian viscosity involving a logarithmic function. However, the original Newtonian constant viscosity can also be replaced by other, even simpler generalized Newtonian models to capture the shear-thinning behavior of blood. Such variable viscosity generalization has been adopted by several authors, including our studies (e.g., [16,17,18]) leading to the generalized Oldroyd-B model with shear-thinning viscosity given by the Cross model with appropriate physiological parameters, as specified below in this paper (Section 3). These generalized Oldroyd-B models however have some limitations. For example, the relaxation times in Oldroyd type models (even generalized) do not depend on the shear rate, which in many cases does not agree with experimental observations and more complex models have been developed for such fluids (see, e.g., [19,20,21] and references cited therein).

The approach used in this paper is different from what is typically adopted when a shear-thinning viscoelastic model is needed. In generalized Oldroyd-B models the standard Oldroyd-B model (predicting constant viscosity) is considered just replacing the model constant viscosity coefficient by a shear rate dependent function (with shear-thinning or shear-thickening behavior).

The procedure we have have used to obtain the shear-thinning behavior in viscoelastic fluids is based on a well known fact that if in the Johnson–Segalman model family the convected derivative is chosen other than upper- or lower-convected, i.e., for the choice of parameter $a\in (-1;1)$, the model will exhibit the shear-thinning behavior, despite keeping constant all the model parameters.

Our aim in this paper is to present a preliminary numerical study of the rate type shear-thinning viscoelastic Johnson–Segalman model in order to better understand its behavior in some practically relevant situation—the blood flow in stenosed vessel. The obtained results are compared to those from generalized Oldroyd-B model, where the shear dependent Cross viscosity function is artificially introduced to capture the shear-thinning characteristic. As already mentioned, this model has been used for blood flow simulations in some of our previous works [16,17,22].

A smooth idealized stenosed blood vessel geometry is used in the numerical simulations, for a finite volume discretization along with an explicit Runge–Kutta time marching scheme, applied to both the Johnson–Segalman and the generalized Oldroyd-B models, at three different flow rates (Reynolds number less than 100). A range of the convected derivative parameter $a\in (-1;1)$ is explored, namely $a=-0.9,-0.8,-0.7$ (leading to different levels of shear-thinning), knowing that the constitutive equation with $a=-1$ corresponds to the Oldroyd-B model with constant viscosity. Numerical results have been obtained to illustrate and provide comparative analysis of the flow dynamics on the generated flow patterns of the pressure, axial velocity and radial velocity for both models in different flow regimes.

2. The Mathematical Model

The equations for the balance of linear momentum and conservation of mass (or incompressibility condition) for isothermal flow are given by:

$$\begin{array}{ccc}\hfill \rho \frac{d\mathit{u}}{dt}& =& \mathrm{div}\mathsf{T}-\nabla p\hfill \end{array}$$

(1)

$$\begin{array}{ccc}\hfill \mathrm{div}\mathit{u}& =& 0\hfill \end{array}$$

(2)

where $\mathit{u}$ is the velocity field, $\rho$ is the constant density of the fluid, p is an isotropic pressure, $d/dt=\partial /\partial t+\mathit{u}·\nabla$ denotes the material time derivative and $\mathsf{T}$ is the extra-stress (or deviatoric stress) tensor that should be defined by a specific constitutive equation relating the state of stress to the kinematic variables, namely the rate of deformation of fluid elements. It accounts for differences in behavior from a purely inviscid incompressible fluid to the simple viscous Newtonian model where $\mathsf{T}=2\mu \mathsf{D}$, $\mu$ being the dynamic viscosity and $\mathsf{D}$ the rate of deformation tensor, defined as the symmetric part of the velocity gradient $\mathsf{D}=(\nabla \mathit{u}+\nabla {\mathit{u}}^T)/2$. Replacing this expression of $\mathsf{D}$ in the previous system, we obtain the classical Navier–Stokes equations for an incompressible viscous fluid.

This work is concerned with the general class of nonlinear rate-type viscoelastic models, of Oldroyd type also called Johnson–Segalman model with a constitutive relation for $\mathsf{T}$ defined by:

$$\mathsf{T}+{\lambda }_1\frac{\delta \mathsf{T}}{\delta t}=2\mu \left(\mathsf{D}+{\lambda }_2\frac{\delta \mathsf{D}}{\delta t}\right)$$

(3)

Here the parameter ${\lambda }_1$ is the relaxation time, the material coefficient ${\lambda }_2$ is the retardation time, with $0\le {\lambda }_2<{\lambda }_1$ and $\mu$ is the dynamical viscosity. The operator $\frac{\delta ·}{\delta t}$ stands for the so-called convected derivative, which is a generalization of the time-derivative chosen so that the principle of frame indifference is verified, meaning that the model is objective under a superposed rigid body motion and the resulting second order tensor is symmetric. Different choices of the convected derivative lead to different models including, e.g., the classical Oldroyd A and B models, as detailed below.

The Johnson–Segalman model can be derived from the simplest rate-type viscoelastic Maxwell model (see, e.g., [23]) with an additional viscosity

$$\mathsf{T}+{\lambda }_1\frac{\delta \mathsf{T}}{\delta t}=2\mu \mathsf{D}$$

(4)

Simple mechanical models can be used to illustrate the typical behavior of viscoelastic materials, where a dashpot (piston moving inside a cylinder filled with a liquid) represents a viscous (Newtonian) fluid and a spring stands for an elastic (Hookean) solid. These elements can be connected in series or in parallel and combined to represent several deformation-stress models to analyze the behavior of different viscoelastic materials [24]. The combination of the Newtonian model and the Maxwell model joined in parallel shown in Figure 1 represents the mechanical analogue of the Johnson–Segalman model (3).

In the Johnson–Segalman model the viscosity $\mu$ is defined as $\mu ={\mu }_s+{\mu }_e$ where ${\mu }_s$ and ${\mu }_e$ are the solvent and the elastic viscosity coefficients. Moreover, parameters ${\lambda }_1$ and ${\lambda }_2$ are such that:

$${\lambda }_1=\frac{{\mu }_e}{G}{\lambda }_2={\lambda }_1\frac{{\mu }_s}{{\mu }_s+{\mu }_e}$$

(5)

where G is the elastic modulus.

In the mechanical analogue scheme shown in Figure 1 we identify two branches. The upper one with dashpot (viscosity ${\mu }_s$) corresponds to the Newtonian solvent fluid behavior and the lower one with dashpot (viscosity ${\mu }_e$) and spring (elastic modulus G) in serial combination, corresponds to the viscoelastic Maxwell fluid. The total stress $\tau ={\tau }_s+{\tau }_e$ is decomposed into the Newtonian solvent contribution ${\tau }_s$ and its viscoelastic complement ${\tau }_e$. Additionally, $\gamma ={\gamma }_s={\gamma }_e$ corresponds to the rate of deformation tensor. Finally, the extra-stress tensor in Equation (3) is decomposed into its Newtonian ${\mathsf{T}}_s$ and elastic ${\mathsf{T}}_e$ parts, where ${\mathsf{T}}_e$ satisfies a constitutive equation of Maxwell type (4) and we get the following relations:

$$\mathsf{T}={\mathsf{T}}_s+{\mathsf{T}}_e$$

(6)

with

$$\begin{array}{ccc}\hfill {\mathsf{T}}_s& =& 2{\mu }_s\mathsf{D}\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {\mathsf{T}}_e+{\lambda }_1\frac{\delta {\mathsf{T}}_e}{\delta t}& =& 2{\mu }_e\mathsf{D}\hfill \end{array}$$

(8)

Convected Derivatives

The general expression of the objective three-parametric family of convected derivatives of any tensor $\mathsf{M}$ is given by:

$${\left(\frac{\delta \mathsf{M}}{\delta t}\right)}_{abc}=\frac{d\mathsf{M}}{dt}-\mathsf{W}\mathsf{M}+\mathsf{M}\mathsf{W}+a(\mathsf{D}\mathsf{M}+\mathsf{M}\mathsf{D})+b(\mathsf{D}:\mathsf{M})\mathsf{I}+c\left(\mathsf{D}\mathrm{tr}\mathsf{M}\right)$$

(9)

where $\mathsf{W}$ represents the anti-symmetric part of the velocity gradient and $\mathsf{I}$ is the identity matrix. In this expression the a, b and c are real parameters, usually with $b=0$ and $c=0$. Different choices of reference frames suitably fixed to the body yield different objective derivatives. Some of the most frequently used derivatives listed in

Table 1. Commonly used convected derivatives.

Name	Notation	Definition	Parameter a
Lower-convected	$\stackrel{△}{\mathsf{M}}$	$\frac{d\mathsf{M}}{dt}+{\mathsf{L}}^T\mathsf{M}+\mathsf{M}\mathsf{L}$	1
Upper-convected	$\stackrel{▿}{\mathsf{M}}$	$\frac{d\mathsf{M}}{dt}-\mathsf{L}\mathsf{M}-\mathsf{M}{\mathsf{L}}^T$	$-1$
Co-rotational (Jaumann)	$\stackrel{\circ }{\mathsf{M}}$	$\frac{d\mathsf{M}}{dt}-\mathsf{W}\mathsf{M}+\mathsf{M}\mathsf{W}$	0

presented below, where $\mathsf{L}=2\mathsf{D}$.

The one-parametric family of convected derivatives can be written as:

$${\left(\frac{\delta \mathsf{M}}{\delta t}\right)}_a=\frac{d\mathsf{M}}{dt}-\mathsf{W}\mathsf{M}+\mathsf{M}\mathsf{W}+a(\mathsf{D}\mathsf{M}+\mathsf{M}\mathsf{D})a\in [-1;1]$$

(10)

This is sometimes named as a Gordon–Schowalter derivative with parameter $a=\xi -1$ where $\xi$ is called slip parameter.

The special cases $a=1$ and $a=-1$ correspond to the classical lower-convected Oldroyd-A and upper-convected Oldroyd-B models, respectively. Usually attention is not focused on the Oldroyd-A models because their viscometric functions do not match the behavior of real fluids.

The Johnson–Segalman model (3) forms a subset of the 8-constant Oldroyd models developed and analyzed in [23]. Roughly speaking this subset could be obtained from the full model (9) by taking the convected derivative in the form (10).

The convected derivative is defined (for any tensor $\mathsf{M}$) by:

$$\frac{\delta \mathsf{M}}{\delta t}=\stackrel{\circ }{\mathsf{M}}+a(\mathsf{D}\mathsf{M}+\mathsf{M}\mathsf{D})a\in [-1;1]$$

(11)

The formula written above (11) is identical to (10) where the definition of the co-rotational derivative $\stackrel{\circ }{\mathsf{M}}$ from Table 1 was used. Using this notation the Johnson–Segalman model [25,26] can be rewritten as follows:

$$\mathsf{T}+{\lambda }_1\stackrel{\circ }{\mathsf{T}}+a{\lambda }_1(\mathsf{D}\mathsf{T}+\mathsf{T}\mathsf{D})=2\mu \left(\mathsf{D}+{\lambda }_2\stackrel{\circ }{\mathsf{D}}+2a{\lambda }_2{\mathsf{D}}^2\right)$$

(12)

3. The Shear-Thinning Viscosity Behavior of the Johnson–Segalman Model

The classical Oldroyd-A and Oldroyd-B models do not predict the viscosity shear-thinning effect. However models based on the Oldroyd-B constitutive equations are quite often used to build generalized models for shear-thinning fluids.The shear-rate dependent apparent viscosity is introduced by generalizing Oldroyd model, where the constant viscosity $\mu$ is replaced by a suitable shear-rate dependent viscosity function $\mu \left(\dot{\gamma }\right)$. This is the classical way of construction of shear-thinning viscoelastic models. The general concept used by most of them is based on viscosity functions of the following form

$$\mu \left(\dot{\gamma }\right)={\mu }_{\infty }+({\mu }_0-{\mu }_{\infty })F\left(\dot{\gamma }\right)$$

(13)

Here ${\mu }_0$ and ${\mu }_{\infty }$ are the asymptotic viscosity values for low and high shear rates. The appropriate transition between these values is carried out by the shear-rate dependent function $F\left(\dot{\gamma }\right)$. This function should satisfy the following natural limit conditions:

$$\underset{\dot{\gamma }\to 0^{+}}{lim}F\left(\dot{\gamma }\right)=1\&\underset{\dot{\gamma }\to \infty }{lim}F\left(\dot{\gamma }\right)=0$$

There are many possible choices for such a function $F\left(\dot{\gamma }\right)$. One of the most frequently used is the generalized Cross formula used in the following model:

$$\mu \left(\dot{\gamma }\right)={\mu }_{\infty }+\frac{{\mu }_0-{\mu }_{\infty }}{{\left(1+{\left(\lambda \dot{\gamma }\right)}^b\right)}^a}$$

(14)

Model adjustable parameters should be obtained by calibration against suitable experimental data set. The following parameters have been used for blood flow simulations based on this model in ([27]):

${\mu }_0=1.6·{10}^{-1}$ Pa· s	${\mu }_{\infty }=3.6·{10}^{-3}$ Pa· s
$a=1.23$, $b=0.64$ s	$\lambda =8.2$ s

This is of course just one of many possible models being used to describe the shear-thinning viscosity of blood. Other choices can for example be found in the book [28], namely in the chapter [29].

Here we propose another approach which allows to take into account the shear-thinning behavior of the fluid, without introducing any artificially generalized viscosity coefficient.

While for $a=-1$ (Oldroyd-B model) and for $a=1$ (Oldroyd-A model) the viscosity is constant, as mentioned above, it is possible to show that any choice of parameter $a\in (-1;1)$ leads to shear-thinning viscoelastic models. This can be simply obtained by reduction of the general expression developed for the 8-constant Oldroyd model in [23]. For the model in the form (12) we get the apparent simple-shear viscosity:

$$\tilde{\mu }\left(\dot{\gamma }\right)=\frac{\mu \left(\dot{\gamma }\right)}{\mu }=\frac{1+\left(\frac{{\lambda }_2}{{\lambda }_1}\right)[1-a^2]{\lambda }_1^2{\dot{\gamma }}^2}{1+[1-a^2]{\lambda }_1^2{\dot{\gamma }}^2}.$$

(15)

In a similar way we obtain the first and second normal stress coefficients ${\mathsf{\Psi }}_1$ and ${\mathsf{\Psi }}_2$

$$\begin{array}{ccc}\hfill {\mathsf{\Psi }}_1& =& 2{\lambda }_1\mu \left[\tilde{\mu }\left(\dot{\gamma }\right)-\left(\frac{{\lambda }_2}{{\lambda }_1}\right)\right]\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\hfill {\mathsf{\Psi }}_2& =& -(1+a){\lambda }_1\mu \left[\tilde{\mu }\left(\dot{\gamma }\right)-\left(\frac{{\lambda }_2}{{\lambda }_1}\right)\right]=-\frac{(1+a)}{2}{\mathsf{\Psi }}_1.\hfill \end{array}$$

(17)

From Formula (15) for the apparent viscosity we conclude that it contains 3 parameters ${\lambda }_1$, a and $\left(\frac{{\lambda }_2}{{\lambda }_1}\right)$. It is easy to see, that for $a^2=1$ the apparent viscosity is constant, while for any other admissible value $a\in (-1;1)$ the $\tilde{\mu }\left(\dot{\gamma }\right)$ will monotonically decrease with growing $\dot{\gamma }$. This behavior can be observed from Figure 2 where the apparent viscosity $\tilde{\mu }\left(\dot{\gamma }\right)$ is defined for ${\lambda }_2/{\lambda }_1=1/100$.

The limit values of $\tilde{\mu }\left(\dot{\gamma }\right)$ are:

$$\underset{\dot{\gamma }\to 0}{lim}\tilde{\mu }\left(\dot{\gamma }\right)=1\mathrm{and}\underset{\dot{\gamma }\to \infty }{lim}\tilde{\mu }\left(\dot{\gamma }\right)=\frac{{\lambda }_2}{{\lambda }_1}$$

(18)

This can be rewritten using (5) and the decomposition $\mu ={\mu }_s+{\mu }_e$, leading to the asymptotic dimensional viscosities ${\mu }_0$ and ${\mu }_{\infty }$ given by:

$${\mu }_0=\mu ={\mu }_s+{\mu }_e\mathrm{and}{\mu }_{\infty }=\mu \frac{{\lambda }_2}{{\lambda }_1}=\mu \frac{{\mu }_s}{{\mu }_s+{\mu }_e}={\mu }_s$$

(19)

From this result it is easy to understand the role of parameter $\alpha =\frac{{\lambda }_2}{{\lambda }_1}$. Setting ${\lambda }_2=0$ we obtain the convected Maxwell model which leads to zero limit viscosity ${\mu }_{\infty }=0$. This is fixed by taking $\alpha \in (0;1)$ which allows to set properly the high shear rate limit viscosity. This behavior is shown in Figure 3, where $\tilde{\mu }\left(\dot{\gamma }\right)$ is plotted for the co-rotational model with $a=0$ (Table 1).

The role of parameters ${\lambda }_1$ and a is slightly more difficult to identify. In Formula (15) these parameters appear together in the term $\left([1-a^2]{\lambda }_1^2\right)$. The appropriate value of this product could be obtained by fitting the expression (15) against some simple-shear viscosity data. The graph of this joint parameter $\left([1-a^2]{\lambda }_1^2\right)$ depending on single arguments a and ${\lambda }_1$ is shown in Figure 4.

Contours of this function are presented in the following Figure 5.

The separate values of ${\lambda }_1$ and a could only be obtained from some additional viscometric data fitting. One possibility is to use the normal stress difference coefficients ${\mathsf{\Psi }}_1$ and ${\mathsf{\Psi }}_2$ to quantify both ${\lambda }_1$ and a. For example, from known ${\mu }_0$ and ${\mu }_{\infty }$ the relaxation time ${\lambda }_1$ can be obtained as:

$$\underset{\dot{\gamma }\to 0}{lim}{\mathsf{\Psi }}_1=2{\lambda }_1\mu \left[1-\left(\frac{{\lambda }_2}{{\lambda }_1}\right)\right]=2{\lambda }_1\left({\mu }_0-{\mu }_{\infty }\right)$$

(20)

In Formula (15) for the apparent simple-shear viscosity, the parameter ${\lambda }_1$ acts as a scaling factor for the argument $\dot{\gamma }$. The effect of varying ${\lambda }_1$ on $\tilde{\mu }\left(\dot{\gamma }\right)$ is shown in Figure 6.

The shear-thinning effect leads to a reduction of the apparent elastic (or polymeric) viscosity. Due to this fact, there is an important difference between the classical characteristic Weissenberg number of the solved problem and the local, effective value for this number. Thus, it seems to be meaningful to generalize and “localize” the definition of the Weissenberg number to account for the shear-thinning behavior. The possible re-definition could be the following:

$$We=\frac{{\lambda }_{1}U}{L}={\lambda }_1\frac{U}{L}\approx We\left(\dot{\gamma }\right)={\lambda }_1\dot{\gamma }$$

(21)

Here the local shear rate is assumed together with the generalized ${\lambda }_1\left(\dot{\gamma }\right)$ taking into account the shear dependency of the elastic component of the viscosity. The apparent viscosity based on Formula (15) could be rewritten in the dimensional form:

$$\mu \left(\dot{\gamma }\right)={\mu }_{\infty }+\frac{{\mu }_0-{\mu }_{\infty }}{1+[1-a^2]{\lambda }_1^2{\dot{\gamma }}^2}$$

(22)

where the second term represents the apparent elastic viscosity ${\mu }_e$, i.e.,

$${\mu }_e\left(\dot{\gamma }\right)=\frac{{\mu }_0-{\mu }_{\infty }}{1+[1-a^2]{\lambda }_1^2{\dot{\gamma }}^2}=\frac{{{\mu }_e}_{{}_0}}{1+[1-a^2]{\lambda }_1^2{\dot{\gamma }}^2}$$

(23)

Using this generalized elastic viscosity and relations (5), the local Weissenberg number can be computed as:

$$We\left(\dot{\gamma }\right)={\lambda }_1\left(\dot{\gamma }\right)\dot{\gamma }=\frac{{\mu }_e\left(\dot{\gamma }\right)}{G}=\frac{{\mu }_0-{\mu }_{\infty }}{G}\frac{\dot{\gamma }}{1+[1-a^2]{\lambda }_1^2{\dot{\gamma }}^2}$$

(24)

This leads to the following generalized local Weissenberg number for the proposed class of shear-thinning models:

$$We\left(\dot{\gamma }\right)=\frac{{\lambda }_1\dot{\gamma }}{1+[1-a^2]{\lambda }_1^2{\dot{\gamma }}^2}=\frac{W{e}_{{}_0}}{1+[1-a^2]W{e}_{{}_0}^2}$$

(25)

In the last expression we have used the low shear limit Weissenberg number definition $W{e}_{{}_0}={\lambda }_1\dot{\gamma }$. The difference between $We\left(\dot{\gamma }\right)$ and $W{e}_{{}_0}$ is that the former one takes into account the shear-thinning, while the later one does not. Due to the shear-thinning viscosity, the local Weissenberg number becomes a nonlinear function of the shear-rate $\dot{\gamma }$. The shear-thinning induced deviation of the Weissenberg number from the classical linear dependency can be observed in Figure 7 and Figure 8 for different values of the convected derivative parameter a.

In terms of the ratio between the “shear-dependent” Weissenberg number $We\left(\dot{\gamma }\right)$ and its constant viscosity counterpart $W{e}_{{}_0}$ it is even more obvious, as observed in Figure 8.

It could easily be shown from (25) that for the Oldroyd-B and Oldroyd-A models, which have no shear-thinning viscosity, the original $We\left(\dot{\gamma }\right)={\lambda }_1\dot{\gamma }$ is recovered. The low-shear limit Weissenberg number leads exactly to the same original linear behavior.

One of the important consequences of the local shear-thinning re-definition of the Weissenberg number is the existence of its local extrema. This is in contrast with the classical, constant viscosity point of view, where the Weissenberg number is a linear (and therefore monotone) function of the shear rate. Based on (25) it could easily be shown that the maximum is attained at:

$${\left[{\lambda }_1\dot{\gamma }\right]}_{max}={\left[W{e}_{{}_0}\right]}_{max}=\frac{1}{\sqrt{1-a^2}}$$

(26)

and consequently the effective Weissenberg number value is given by

$$We\left({\left[{\lambda }_1\dot{\gamma }\right]}_{max}\right)=\frac{1}{2}\frac{1}{\sqrt{1-a^2}}=\frac{1}{2}{\left[{\lambda }_1\dot{\gamma }\right]}_{max}$$

(27)

4. Numerical Tests

The numerical solution of the above described model is based on a finite-volume semi-discretization in space on structured grids and an explicit Runge–Kutta time integration scheme. Steady solution is searched by a time-marching approach, i.e., the unsteady governing system is solved with steady boundary conditions and stationary solution is recovered when $t⟶\infty$. Details of the numerical method can be found in our previous works [16,17,18].

4.1. Computational Geometry and Mesh

In this study, we follow the experimental set-up introduced in [27]. The first computational domain represents a non-symmetric cosine-shaped vessel stenosis shown in Figure 9.

The vessel is three-dimensional, rotationally symmetric with maximum diameter $D=6.2$ mm which reduces to its one half in the stenosed region. This leads to a 4:1 cross-sectional area reduction and thus to a significant local flow acceleration. The effect of viscoelasticity in this region will be studied in detail.

We have used structured, wall fitted mesh with hexahedral cells. Multiblock mesh structure was adopted to avoid high distortion of cells. Figure 10 shows the actual grid used in all presented simulations.

The grid had $16\times 48+12\times 12=912$ cells in each crossection ($z=const$), where the smallest (near wall) cell size was set to $1\%$ of the vessel radius, i.e., $\Delta x=0.01D/2$. From this smallest scale the cells are gradually (an smoothly) increased towards the inner grid block surrounding the channel axis. In the axial direction 143 cells were used, from which 80 equispaced cells were placed in the channel contraction. This fine cells in the stenosis are again gradually enlarged towards the channel inlet and outlet. This grid provided consistent results for all presented cases and further refinement did not led to any visible changes with respect to presented simulations.

4.2. Model Parameters

4.2.1. The Generalized Oldroyd-B Model

The basic set of parameters for the generalized (shear-thinning) Oldroyd-B model was taken from [27] and previously used in some of our studies [16,17,18]. The shear-thinning viscosity is given by the Cross model (14):

$$\mu \left(\dot{\gamma }\right)={\mu }_{\infty }+\frac{{\mu }_0-{\mu }_{\infty }}{{\left(1+{\left(\lambda \dot{\gamma }\right)}^b\right)}^a}$$

with the following parameters suitable for blood flow simulations [27]:

${\mu }_0=1.6·{10}^{-1}$Pa· s	${\mu }_{\infty }=3.6·{10}^{-3}$ Pa· s
$a=1.23$, $b=0.64$ s	$\lambda =8.2$ s

The data for the viscoelastic part of the model are taken from the same paper [27], where:

${\mu }_e=4.0·{10}^{-4}$Pa· s	${\mu }_s=3.6·{10}^{-3}$Pa· s
${\lambda }_1=0.06$ s	${\lambda }_2=0.054$ s
${\mu }_0=\mu ={\mu }_s+{\mu }_e$	$\rho =1050$ kg· m${}^{-3}$

This generalized Oldroyd-B model is hereafter denoted as GOB model.

4.2.2. The Johnson–Segalman Model

Based on these parameters, the shear-thinning Johnson–Segalman model was adjusted. Namely the viscosities ${\mu }_0$ and ${\mu }_{\infty }$ were preserved together with the retardation time scale ${\lambda }_2$. Using the viscosity relations (19) we get the appropriate solvent and elastic viscosity ${\mu }_s$ and ${\mu }_e$, respectively

$${\mu }_s=\mu \frac{{\lambda }_2}{{\lambda }_1}=\mu \frac{{\mu }_s}{{\mu }_s+{\mu }_e}={\mu }_{\infty };{\mu }_e=\mu -{\mu }_s={\mu }_0-{\mu }_{\infty }$$

(28)

leading to the relaxation time ${\lambda }_1$ given by:

$${\mu }_s=\mu \frac{{\lambda }_2}{{\lambda }_1}⟹\frac{{\mu }_s}{\mu }=\frac{{\mu }_{\infty }}{{\mu }_0}=\frac{{\lambda }_2}{{\lambda }_1}⟹{\lambda }_1={\lambda }_2\frac{{\mu }_0}{{\mu }_{\infty }}$$

(29)

This setting will guarantee that the asymptotic shear viscosity of the Johnson–Segalman model is the same as for the standard GOB model. The only remaining parameter $a\in (-1;1)$ in the Gordon–Schowalter convected derivative will be independently adjusted to assess its effect on the model predictions. Values will be chosen close to the Oldroyd-B model (upper convected derivative), which means close to $a=-1$.

5. Numerical Results

In this paper, only some initial simulation results are presented to demonstrate the qualitative comparison and agreement of the widely used shear-thinning generalized Oldroyd-B model and the proposed Johnson–Segalman model in the shear-thinning regime with $\left|a\right|<1$. The standard GOB model is used to obtain reference flow fields in the axisymmetric stenosis for three different flow rates $Q=0.5,1.0,2.0$ cm${}^3$/s. The resulting fields of axial velocity, radial velocity and pressure are compared (for the same flow rates) with those obtained using Johnson–Segalman model with a set of parameters a equal to $a=-0.9,-0.8,-0.7$.

5.1. Generalized Oldroyd-B Model—Reference Results

The flow in the axisymmetric stenosis with circular cross-section is simulated for three different flow regimes. In this section, the simulation results are grouped together according to the flow rate. The contour fields are shown in the longitudinal axial cut, showing always the pressure, axial velocity and radial velocity fields, including the corresponding color scales. For each one of the flow rates different color scales were chosen to better capture the individual flow features. Scales were kept dimensional, since simple non-dimensional scaling based, e.g., on flow rates are not optimal due to strongly nonlinear viscosity behavior.

In all presented simulations the pressure is set to zero at the outlet (right) boundary, therefore the maximum is reached at the inlet (left) boundary. The density of pressure contour lines in the stenosed region indicates high local pressure gradient leading to significant flow acceleration visible in the axial velocity contours (central part of the stenosed region). The radial velocity field is axially anti-symmetric, showing the appropriate flow tendency towards the axis or the wall, following the same behavior of the nearest wall geometry.

Comparing the profiles of the velocity fields in Figure 11 for $Q=0.5$ cm${}^3$/s with results obtained for higher flow rates shown in Figure 12 and Figure 13 reveals that for lower flow rates the velocity profiles become smoother, with more rounded contour lines downstream the stenosis, due to lower flow separation and recirculation. This is caused by the effect of shear-thinning behavior, where for lower shear rates (i.e., lower flow rates) the shear dependent viscosity provides higher values, closer to ${\mu }_0$. In a straight tube (non-stenosed vessel) the zero shear rate is always achieved at the central axis and its value increases towards the wall (see, e.g., [18,22,30]). The viscosity thus behaves accordingly, exhibiting its maximum at the axis, dropping significantly towards the wall. Of course, as demonstrated in [22], the behavior is more complex in non-trivial geometries with significant flow acceleration/deceleration or streamlines curvature.

The profiles of the velocity fields in Figure 11, Figure 12 and Figure 13, have shown the expected behavior, including their steadiness and axial symmetry. Of course the behavior of many other mechanical quantities, such as viscosity, shear rate, stress tensor components, pressure drop or wall shear stress could be represented for these test cases. However, the shown fields are quite representative, in view of the basic qualitative comparison designed for this paper.

5.2. Johnson–Segalman Model—Comparison

The same setup as used for the GOB model, with parabolic velocity profile at the inlet (left) boundary and fixed (zero) pressure at the outlet (right) boundary was also used in the numerical simulations with the Johnson–Segalman model. In this section, the contour fields are primarily grouped according to flow variables and then by flow rates.

The main goal is to compare the results obtained using the Johnson–Segalman shear-thinning model with those obtained using the standard reference GOB model. Therefore in the figures presented below we always group the results from the GOB model with the corresponding predictions of the Johnson–Segalman model obtained by setting the parameter a to $a=-0.9,-0.8,-0.7$.

5.2.1. Axial Velocity

For the considered stenosed vessel test case, the axial velocity is the dominant and most interesting flow quantity. The velocity fields can be directly compared to identify some of the most relevant differences between models predictions. The inlet velocity profile is exactly the same for all models. It is the standard (second order) parabolic profile with maximum velocity in the center, respecting the no-slip velocity condition on the wall. It corresponds to the fully developed velocity profile of a constant viscosity Poiseuille flow of both Newtonian and Oldroyd-B (non-generalized) model. For shear thinning models the velocity profile becomes more flat in the center with steep decay just very close to the wall. Such behavior can be observed when comparing the inlet and outlet profiles in Figure 14, where the outlet velocity profiles are evidently more flat (due to shear-thinning) in the center compared to constant viscosity profiles prescribed at the inlet.

This clearly shows the shear-thinning behavior of the Johnson–Segalman model for $\left|a\right|<1$. For the smallest flow rate $Q=0.5$ cm${}^3$/s shown in Figure 14 the velocity fields downstream the stenosis differ substantially from the GOB model. For the value $a=-0.7$ the Johnson–Segalman model seems to provide the most similar results compared to reference GOB solution, but it’s still quite far. The viscoelastic extra stress probably plays a more important role at this flow rate.

The situation is visibly improved at higher flow rates as it can be observed in Figure 15 for $Q=1.0$ cm${}^3$/s and Figure 16 for $Q=2.0$ cm${}^3$/s. At these flow rates the differences between the reference GOB model solution and the Johnson–Segalman model solution can be minimized by choosing a suitable value of parameter a. The optimal value seems to be close to $a=-0.8$ for $Q=1.0$ cm${}^3$/s (see Figure 15) and closer to $a=-0.9$ for $Q=2.0$ cm${}^3$/s (see Figure 16).

Overall, the axial velocity field comparison shows the expected shear-thinning behavior. However, the model regulation just by the single parameter a seems to be insufficient and optimal values may be dependent on the shear rate (flow rate).

For the generalized Oldroyd-B model (GOB), at higher flow rates ($Q=1.0$ cm${}^3$/s and $Q=2.0$ cm${}^3$/s) a rather large flow separation and recirculation region appears in the expanding part of the stenosis. The zone of reversal flow (where axial velocity is negative) is marked by purple color in Figure 17, showing the extent of the reversal flow predicted by considered models at two higher flow rates (color scale is omitted as it is identical to Figure 15 and Figure 16).

Evidently the setting of parameter a has significant impact on the size of the recirculation zone. This should be considered when looking for the optimal value of a in specific fluid simulations. It should be kept in mind that in all simulations the asymptotic viscosities ${\mu }_0$ and ${\mu }_{\infty }$ are set identically, so the different extent of the recirculation zone cannot easily be explained using the apparent viscosity. Here probably the differences in first and second normal stress differences ${\mathsf{\Psi }}_1$, ${\mathsf{\Psi }}_2$ can come to play as they depend on the parameter a (see (16) and (17)).

5.2.2. Radial Velocity

Concerning the radial velocity, the model change effects are more delicate and less important. For the lowest flow rate $Q=0.5$ cm${}^3$/s the radial velocity fields comparison shown in Figure 18 confirms significant differences of the Johnson–Segalman model predictions with respect to the GOB model reference results. The major ones appear downstream the stenosis near the wall, with a tendency to vanish far from the flow separation and recirculation away from the stenosed region.

As for the axial velocity, the radial velocity fields also become nearly analogous to the reference GOB model solutions obtained at higher flow rates. The Figure 19 and Figure 20 indicate that a suitable choice of the parameter a in the Johnson–Segalman model can lead to an almost identical velocity field as the one obtained for the reference GOB model. Not only the main axial flow patterns are practically equal, but also the secondary flows represented by the radial velocity are very similar.

5.2.3. Pressure Field

A direct comparison of the pressure fields between different models is more complex and less useful. In viscoelastic fluid models, the role and interpretation of the quantity denoted by p may be different compared to what we are used to in the Newtonian models. Physically, pressure is defined as the spherical part of the complete stress tensor, so pressure is proportional to the trace of the full stress tensor and appears as the Lagrange multiplier arising from the incompressibility condition. In our (and in many other) viscoelastic models, part of the physical pressure is contained in the variable (scalar field) p present in the gradient term in the momentum equations. Another part of the physical pressure may be hidden in the trace of the (viscoelastic) stress tensor, which is in general not traceless, in contrast to the Newtonian viscous stress tensor. Therefore, the direct comparison of the pressure fields p does not give the complete information about the physical pressure that can be measured in the flow by appropriate investigations.

Thus, the pressure fields shown below in Figure 21, Figure 22 and Figure 23, correspond just to the part denoted by p in the momentum equations and therefore should be used and interpreted with caution. All we can observe here from the pressure fields is that the general character does not differ substantially from what we have seen for the reference GOB model solution. Keeping the same pressure values at the outlet for all models, we can again observe the maximum values being always achieved at the inlet. The pressure drop (pressure difference between inlet and outlet) visibly changes, however again direct comparison is not possible, because different part of the physical pressure is still hidden in the viscoelastic extra stress tensor.

In any case it is obvious that the pressure fields p for the Johnson–Segalman and the reference GOB solutions become closer at higher flow rates. This can again be attributed to the shear-dependent response of the Johnson–Segalman model leading to extra stress reduction due to shear -thinning. This is confirmed when comparing the results for given flow rate depending on the choice of parameter a.

6. Conclusions and Remarks

Although the Johnson–Segalman model is known for a long time, it’s very seldom used in simulations of shear thinning viscoelastic fluids of practical interest. The overview of basic properties we provide in this paper, together with some initial numerical simulations are aimed at giving a deeper insight and motivation for a broader use of this simple viscoelastic shear-thinning model. The case of blood flow simulation was chosen from many possible examples of suitable applications of the Johnson–Segalman model in the shear-thinning regime. This choice was mainly motivated by our previous experience with other viscoelastic models for this kind of applications [16,22] or [18].

In this paper, we only show a small set of initial results, so we are still quite far from fully understanding the model and its behavior. Nevertheless, already these initial simulations have brought some important information as well as some open questions. The main points are summarized below.

• The Johnson–Segalman model, using the Gordon–Schowalter family of convected derivatives (with $\left|a\right|\ne 1$) naturally exhibits the shear-thinning behavior similar to what is artificially introduced into the generalized Oldroyd-B model. This behavior was documented analytically and confirmed in numerical simulations.
• It was shown how the Johnson–Segalman model parameters can be set and adjusted to mimic the asymptotic behavior of the shear-thinning generalized Oldroyd-B model (i.e., the viscosities ${\mu }_0$, ${\mu }_{\infty }$). The shape of the viscosity-shear $\mu \left(\dot{\gamma }\right)$ dependence for the Johnson–Segalman model is given by (15) and (22). This means the model can only be fine-tuned by choosing its few parameters, but its general behavior cannot be dramatically altered. This might be seen as a disadvantage when compared to the more standard generalized Oldroyd-B models, where the type of the viscosity-shear rate dependence can be freely chosen and adjusted (even for shear-thickening fluids).
• The Johnson–Segalman model predictions seem to follow (at least at higher flow rates) the reference solution obtained from the generalized Oldroyd-B model. It should be kept in mind that numerical predictions of the two models are compared. Therefore, from this purely numerical comparison we cannot assess which of the models is better suited for blood (or any other fluid) flow simulations. For such decision the model predictions must be compared to real fluid experimental data and the level of agreement may differ from one fluid (or solved case) to another.
• The Johnson–Segalman model was used in this paper “as it is”, with no sophisticated fitting of the shear-thinning part. Only the asymptotic viscosities were preserved for both the GOB and Johnson–Segalman models as the key measurable physical parameters of practical interest. Better agreement of results can probably be achieved by fitting the complete viscosity curves of both models. Then some other, more optimal choice of relaxation/retardation time ${\lambda }_1$, ${\lambda }_2$ and parameter a can lead to better mutual agreement of the obtained results.
• Open remains the question of possible physical interpretation of the variation of the parameter a in the convected derivative. As it was shown for $\left|a\right|<1$ it is responsible for the apparent viscosity shear-thinning behavior. However, it’s also worth noting that the Johnson–Segalman model (for $a\ne -1$) predicts non-zero second normal stress difference ${\mathsf{\Psi }}_2\ne 0$ (see (17), which is in contrast with the classical Oldroyd-B model and all its trivial generalizations where ${\mathsf{\Psi }}_2=0$. It appears that the Johnson–Segalman model might be a possible option in cases where the non-zero second normal stress difference was experimentally observed.