Shear rate projection schemes for non-Newtonian fluids

TL;DR

This paper introduces shear rate projection schemes for simulating non-Newtonian fluids, extending existing methods to ensure consistent pressure calculations and demonstrating their accuracy through numerical experiments.

Contribution

It proposes a new family of shear rate projection methods for generalized Newtonian fluids, improving pressure consistency in numerical simulations.

Findings

Methods produce consistent pressure for non-Newtonian fluids.

Numerical experiments validate accuracy with a manufactured solution.

Flow past a cylinder with Carreau model demonstrates practical applicability.

Abstract

The operator splitting approach applied to the Navier-Stokes equations, gave rise to various numerical methods for the simulations of the dynamics of fluids. The separate work of Chorin and Temam on this subject gave birth to the so-called projection methods. The basic projection schemes, either the incremental or non-incremental variant (see [1]) induces an artificial Neumann boundary condition on the pressure. By getting rid of this boundary condition on the pressure, the so-call rotational incremental pressure-correction scheme as proposed by Timmermans et al. [2] for Newtonian fluids with constant viscosity gives a consistent equation for the pressure. In this work we propose a family of projection methods for generalized Newtonian fluids based on an extension of the rotational projection scheme. Called shear rate projections, these methods produces consistent pressure when applied to generalized Newtonian fluids. Accuracy of the methods will be illustrated using a manufactured solution. Numerical experiments for the flow past a cylinder, with a Carreau rheological model, will also be presented.

Figures

Tables

Table 1. Table 1: Drag coefficient ( C D subscript 𝐶 𝐷 C_{D} ) and wake lenght ( L 𝐿 L ) at Reynolds R ​ e = 10 𝑅 𝑒 10 Re=10 , Carreau number C U = 10 , 20 subscript 𝐶 𝑈 10 20 C_{U}=10,20 and for a power index ( m 𝑚 m ) ranging from 0.4 0.4 0.4 to 1.0 1.0 1.0 . Results for the steady states mixed (denoted SS-Mix) and unsteady SRP im schemes.

			$C_D$				$L$
$C_U$	$m$		SS-Mix	SRPim	${\Delta }_{C_D}(\%)$		SS-Mix	SRPim	${\Delta }_L(\%)$
10	1.0		2.7517	2.7537	0.07		0.2364	0.2376	0.52
	0.9		2.4697	2.4721	0.10		0.3177	0.3205	0.87
	0.8		2.1955	2.1963	0.04		0.4280	0.4298	0.42
	0.7		1.9292	1.9289	0.01		0.5808	0.5816	0.12
	0.6		1.6747	1.6735	0.07		0.7955	0.7946	0.11
	0.5		1.4348	1.4334	0.10		1.1055	1.1031	0.22
	0.4		1.2140	1.2150	0.08		1.5617	1.5642	0.16
20	1.0		2.7517	2.7537	0.07		0.2364	0.2376	0.52
	0.9		2.3920	2.3908	0.05		0.3703	0.3710	0.19
	0.8		2.0524	2.0528	0.02		0.5627	0.5639	0.22
	0.7		1.7395	1.7404	0.05		0.8437	0.8447	0.11
	0.6		1.4583	1.4588	0.04		1.2652	1.2662	0.08
	0.5		1.2123	1.2118	0.04		1.9086	1.9077	0.05
	0.4		1.0042	1.0045	0.03		2.8702	2.8709	0.02

Equations

$$\sigma =-pI+\nu (\Vert \mathbf{D}(\mathbf{u})\Vert )\mathbf{D}(\mathbf{u}).$$

$$\sigma (\mathbf{u})=-pI+2\nu \mathbf{D}(\mathbf{u}).$$

$$\{\begin{array}{rl}\rho & \frac{\partial \mathbf{u}}{\partial t}+\rho \left(\mathbf{u}\cdot \nabla \right)\mathbf{u}-\nabla \cdot (2\nu \mathbf{D}(\mathbf{u}))+\nabla p=\mathbf{f},\\ & \nabla \cdot \mathbf{u}=0\end{array}\text{mbox}on{\Omega }_t,$$

$$\mathbf{D}(\mathbf{u})=\frac{1}{2}(\nabla \mathbf{u}+\nabla {\mathbf{u}}^t).$$

$$\nu (t,\mathbf{u})=\nu (t,\mathbf{x},\mathbf{u},\mathbf{D}(\mathbf{u}),...)$$

$$\mathbf{u}(0,\mathbf{x})={\mathbf{u}}_0(\mathbf{x})\in L^2(\Omega {)}^d\text{mbox}with\nabla \cdot {\mathbf{u}}_0=0$$

$$\mathbf{u}=0\text{mbox}on{\Gamma }_D\subseteq \partial \Omega .$$

$${\mathcal{W}}_r=\left\{\mathbf{v}\in L^r(0,T;V_r):{\partial }_t\mathbf{v}\in L^{r/(r-1)}(0,T;V_r^{\prime })\right\},\mathcal{M}=L^2(0,T;L^2(\Omega ))$$

$$V_r=\left\{\mathbf{v}\in {\left(W^{1,r}(\Omega )\right)}^d:\nabla \cdot \mathbf{v}=0\text{ on }\Omega ,\mathbf{v}=0\text{mbox}on{\Gamma }_D\right\}.$$

$$\begin{array}[]{c}\nu:(0,T)\times\Omega\times\mathbb{R}^{d}\times\mathbb{R}^{d\times d}\longrightarrow\mathbb{R}\vspace*{4pt}\\ (t,{\boldsymbol{x}},{\boldsymbol{u}},{\boldsymbol{F}})\longmapsto\nu(t,{\boldsymbol{x}},{\boldsymbol{u}},{\boldsymbol{F}})\\ \end{array}$$

$$\nu (t,\mathbf{x},\mathbf{u},\mathbf{F})\le \Vert \mathbf{F}{\Vert }^{-1}\left({\lambda }_1(t,x)+c_1\left(\Vert \mathbf{u}{\Vert }^{r-1}+\Vert \mathbf{F}{\Vert }^{r-1}\right)\right)$$

$$\nu (t,\mathbf{x},\mathbf{u},\mathbf{F})\ge \Vert \mathbf{F}{\Vert }^{-2}\left(-{\lambda }_2(t,x)-{\lambda }_3(t,x)\Vert \mathbf{u}{\Vert }^{\alpha }+c_2\Vert \mathbf{F}{\Vert }^r\right)$$

$$\left(\nu (t,\mathbf{x},\mathbf{u},\mathbf{F})\mathbf{F}-\nu (t,\mathbf{x},\mathbf{u},\mathbf{G})\mathbf{G}\right):(\mathbf{F}-\mathbf{G})\ge 0\forall \mathbf{F},\mathbf{G}\in {\mathbb{R}}^{d\times d}.$$

$$\nu (\mathbf{D}(\mathbf{u}))={\nu }_{\infty }+({\nu }_0-{\nu }_{\infty }){\left(C_0+{\lambda }^2\Vert \mathbf{D}(\mathbf{u}){\Vert }^2\right)}^{\frac{m-1}{2}}m,C_0,{\nu }_{\infty },{\nu }_0\ge 0.$$

$$D_t{\mathbf{u}}^{n+1}=\frac{3\mathbf{u}(t^{n+1},\mathbf{x})-4\mathbf{u}(t^n,\mathbf{x})+\mathbf{u}(t^{n-1},\mathbf{x})}{2\Delta t}.$$

$$\{\begin{array}{rl}\rho & D_t{\mathbf{u}}^{n+1}+\rho \left({\mathbf{u}}^{n+1}\cdot \nabla \right){\mathbf{u}}^{n+1}-\nabla \cdot (2{\nu }^{n+1}({\mathbf{u}}^{n+1})\mathbf{D}({\mathbf{u}}^{n+1}))+\nabla p^{n+1}={\mathbf{f}}^{n+1},\\ & \nabla \cdot {\mathbf{u}}^{n+1}=0\end{array}\text{mbox}on\Omega ,$$

$$\{\begin{array}{rl}\rho & D_t^{\star }\tilde{\mathbf{u}}+\rho \left(\tilde{\mathbf{u}}\cdot \nabla \right)\tilde{\mathbf{u}}-\nabla \cdot (2{\nu }^{n+1}(\tilde{\mathbf{u}})\mathbf{D}(\tilde{\mathbf{u}}))+\nabla p^n={\mathbf{f}}^{n+1}\text{vspace}\ast 4pt\\ & \tilde{\mathbf{u}}=0\text{mbox}on{\Gamma }_D\end{array}$$

$$\{\begin{array}{rl}{\mathbf{u}}^{n+1}& =\tilde{\mathbf{u}}-\frac{2\Delta t}{3\rho }\nabla \phi \\ \nabla \cdot {\mathbf{u}}^{n+1}& =0\\ {\mathbf{u}}^{n+1}\cdot \mathbf{n}& =0\text{mbox}on{\Gamma }_D\end{array}$$

$$p^{n+1}=p^n+\phi$$

$$D_t^{\star }\tilde{\mathbf{u}}=\frac{3\tilde{\mathbf{u}}-4{\mathbf{u}}^n+{\mathbf{u}}^{n-1}}{2\Delta t}=D_t{\mathbf{u}}^{n+1}-\frac{1}{\rho }\nabla \phi .$$

$$\{\begin{array}{rl}\rho & D_t^{\star }\tilde{\mathbf{u}}+\rho \left(\tilde{\mathbf{u}}\cdot \nabla \right)\tilde{\mathbf{u}}-\nabla \cdot (2{\nu }^{n+1}(\tilde{\mathbf{u}})\mathbf{D}(\tilde{\mathbf{u}}))+\nabla p^n={\mathbf{f}}^{n+1}\\ & \tilde{\mathbf{u}}=0\text{mbox}on{\Gamma }_D\end{array}$$

$$\{\begin{array}{rl}{\mathbf{u}}^{n+1}& =\tilde{\mathbf{u}}-\frac{2\Delta t}{3\rho }\nabla \phi \\ \nabla \cdot {\mathbf{u}}^{n+1}& =0\\ {\mathbf{u}}^{n+1}\cdot \mathbf{n}& =0\text{mbox}on{\Gamma }_D\end{array}$$

$$p^{n+1}=p^n+\phi +\psi$$

$$\nabla p^n=\nabla p^{n+1}-3\rho \left(\frac{\tilde{\mathbf{u}}-{\mathbf{u}}^{n+1}}{2\Delta t}\right)-\nabla \psi .$$

$$\rho D_t{\mathbf{u}}^{n+1}$$

$$\nabla \psi =\nabla \cdot (2{\nu }^{n+1}({\mathbf{u}}^{n+1})\mathbf{D}({\mathbf{u}}^{n+1}))-\nabla \cdot (2{\nu }^{n+1}(\tilde{\mathbf{u}})\mathbf{D}(\tilde{\mathbf{u}})).$$

$$\Delta \psi =\nabla \cdot \nabla \cdot \left(2({\nu }^{n+1}({\mathbf{u}}^{n+1})-{\nu }^{n+1}(\tilde{\mathbf{u}}))\mathbf{D}(\tilde{\mathbf{u}})\right)-\nabla \cdot \nabla \cdot \left(2{\nu }^{n+1}({\mathbf{u}}^{n+1})\mathbf{D}(\tilde{\mathbf{u}}-{\mathbf{u}}^{n+1})\right).$$

$$\nabla p^{n+1}\cdot \mathbf{n}=\left({\mathbf{f}}^{n+1}-\rho D_t{\mathbf{u}}^{n+1}+\rho \left(\tilde{\mathbf{u}}\cdot \nabla \right)\tilde{\mathbf{u}}-\nabla \cdot (2{\nu }^{n+1}({\mathbf{u}}^{n+1})\mathbf{D}({\mathbf{u}}^{n+1}))\right)\cdot \mathbf{n}.$$

$$\{\begin{array}{rl}& 3\rho \frac{\tilde{\mathbf{u}}}{2\Delta t}+\rho \left({\tilde{\mathbf{u}}}^{\star }\cdot \nabla \right)\tilde{\mathbf{u}}-\nabla \cdot (2{\nu }^{\star }\mathbf{D}(\tilde{\mathbf{u}}))={\mathbf{f}}^{n+1}-\nabla p^n+\rho \frac{4{\mathbf{u}}^n-{\mathbf{u}}^{n-1}}{2\Delta t}\\ & \tilde{\mathbf{u}}=\text{0}\text{mbox}on{\Gamma }_D\end{array}$$

$$\{\begin{array}{rl}& 3\rho \frac{{\tilde{\mathbf{u}}}_{k+1}}{2\Delta t}+\rho \left({\tilde{\mathbf{u}}}^{\star }\cdot \nabla \right){\tilde{\mathbf{u}}}_{k+1}-\nabla \cdot (2{\nu }^{\star }\mathbf{D}({\tilde{\mathbf{u}}}_{k+1}))={\mathbf{f}}^{n+1}-\nabla p^n+\rho \frac{4{\mathbf{u}}^n-{\mathbf{u}}^{n-1}}{2\Delta t}\\ & {\tilde{\mathbf{u}}}_{k+1}=\text{0}\text{mbox}on{\Gamma }_D.\end{array}k=1,2,\dots$$

Full text

††journal:

1 Introduction

Solving the Navier–Stokes system describing unsteady flows is a theoretical and numerical challenge for a homogeneous fluid and, a fortiori, for heterogeneous fluids. To numerically simulate such flows, apart from the heterogeneous nature of the fluid, we are confronted with two major difficulties : very complex dynamics, needing accurate time approximation (high order time approximation or fine time step); complex spatial behaviour that dictates the use of very fine meshes when solving the incompressible Navier–Stokes equations. The construction of efficient solver to achieve good approximation in a reasonable computational time is a difficult task. Since we are interested in three dimensional problems, numerical methods based on mixed formulation (for example [3, 4] in the finite element case), although precise, can be time and resources consuming.

In designing an effective time-marching techniques for this problem we need to address the fact that the incompressibility constraint in the Navier–Stokes system gives the problem a saddle point structure. One way to overcome this structure is found in the pioneering works of Chorin [5, 6] and Temam [7] who introduced projection methods. The idea is to apply a fractional time step to decouple the incompressibility constraint from the diffusion operator based on the Helmoltz decomposition (see [4] for instance).

The projection method as originally presented has some drawbacks: limited precision of the resulting algorithm (see Rannacher [8] and Shen [9], the use of any higher-order time stepping scheme does not improve the overall accuracy), and artificial boundary condition of the resulting pressure.

To alleviate those deficiencies, numerous variants have been proposed over the years ([10, 8, 11, 12, 13]). For an interesting overview, we also refer to Guermond and coauthors [1, 13]. Notably, those limitations lead to the introduction of an incremental projection scheme (proposed by Goda et al [14, 15]), a projection where the viscous equation takes into account the pressure at the previous time-step. In the context of homogeneous fluids, this variation of the original scheme is commonly used since it provides an improved precision (see [1]).

Timmermans et al. [2] proposed the rotational projection scheme (RP) or rotational incremental pressure-correction scheme in Guermond’s terminology. This scheme gives a consistent boundary condition for the pressure and better precision for the velocity and pressure, see for instance Guermond et al. [1]. However, this approach is only valid for homogeneous viscosity as illustrated in [16]; limiting its use to homogeneous Newtonian flows. Still based on the original idea of Chorin and Temam, the more general shear rate projection (SRP) proposed by Deteix et al. [16] allows the treatment of heterogeneous viscosity and improves the accuracy of the incremental projection (the rate of convergence are comparable to those of the rotational projection for Newtonian fluids). This enhanced precision makes the shear rate projection attractive in applications such as natural convection, Allen-Cahn or Cahn-Hilliard flows but also non–Newtonian fluids.

Non–Newtonian model can be summarily described as modelling complex flows using constitutive models which involve dependency of the local stress on the velocity gradient (shear-rate dependant fluid) and possibly the deformation history of the fluid. In this work we are interested in generalized Newtonian fluids. Generalized Newtonian fluids form a subclass of non–Newtonian models. They are the simplest extension from Newtonian models to non–Newtonian models. They model phenomenon where the flow affect the viscosity of the fluid, but does not change the ”Newtonian nature” of the constitutive law :

[TABLE]

The velocity fields of a generalized Newtonian fluids is instantaneously modified by the stresses (the history of the flow has no effect). These fluids are characterized by the derivative of the viscosity with respect to $\|{\bf D}(u)\|$ (the tensorial norm). Monotone viscosity leads to two possible states: shear-thickening where resistance to shear increases as the shear rate increases and shear-thinning fluids having the opposite behaviour. Of course more complex behaviour can be modelled with non-monotone relation between shear-rate and viscosity.

The constitutive law is purely phenomenological as it try to replicate the stresses due to the applied flow by a shear-dependent effective viscosity. This give rise to various phenomenological description of viscosity adapted to various empirical observations. Carreau–Yasuda, power law, Cross, etc. [17, 18] are constitutive equations for $\nu$ of particular interest as they seems to be frequently used in numerous situations. There is an extensive and recent literature ([19, 20, 21, 22, 23] for example) concerning existence and regularity of solutions of the generalized Newtonian flow problem for multiple constitutive laws.

Our goal is to propose an original projection method for the numerical simulation of generalized Newtonian flow. To achieve this, we will restrict ourselves to a rheological and theoretical context insuring validity of the continuous model (more precisely we refer to [21]). The shear rate projection as proposed in [16] does not take into account the explicit dependence of the viscosity upon velocity. However, the same paradigm applies to the projection method when used on generalized Newtonian fluid: some information, easily extracted from the predicted velocity, is ignored. We propose a modification of the shear rate projection, recuperating part of the lost information, making it possible to gain in accuracy and consistency. Based on this generalized version of the SRP, a family of numerical scheme based on the finite element method combined to a second order time approximation is proposed.

The last section will be devoted to explore different representative of this family of methods. The implementation is based on the finite element library FreeFem, [24], code efficiency and optimal strategies will be minimal and performance measure will be limited to comparison of physical quantities for various methods. Time and spatial accuracy will be illustrated first, followed by a classical application, that can be regarded as benchmark: the steady flow past a cylinder.

2 Problem setting

We consider an unsteady flow of a generalized Newtonian fluid : a fluid having a non-homogeneous viscosity $\nu$ depending on time, space, velocity or shear rate of the fluid, possibly the fluid pressure, temperature or other external quantities and a constitutive equations of the form

[TABLE]

Assuming a Boussinesq-type model, the density variations are neglected with the exception of the forces term. These assumptions leads to the Navier–Stokes equations for the fluid velocity ${\bf u}$ and pressure $p$

[TABLE]

where the shear rate is defined as

[TABLE]

The general expression

[TABLE]

represents the viscosity of the fluid and ${\boldsymbol{f}}$ represents external volumic forces (such as gravity). System (1) is completed with the following initial data:

[TABLE]

and we consider homogeneous Dirichlet boundary conditions on a non-empty, but possibly limited, part of the boundary of the domain

[TABLE]

2.1 Existence and regularity of solutions

Let $\Omega$ be a smooth domain in $\mathbb{R}^{d},~{}d=$ 2 or 3 which satisfy the inf–sup conditions (see [4]). Let $\partial\Omega=\Gamma_{D}\cup\Gamma_{N}$ the boundary of $\Omega$ (possibly $\Gamma_{N}=\emptyset$)), ${\boldsymbol{n}}$ is the exterior normal vector and $\Omega_{t}$ the open set $\Omega\times\left(0,T\right)$, where $T>0$ is the final time.

We denote by $L^{2}(\Omega)$ the space of square integrable functions defined on $\Omega$. The Sobolev spaces $W^{p,q}(\Omega)$, $q\geq 0$, are the spaces of functions in $L^{q}(\Omega)$ with generalized partial derivatives belonging to $L^{q}(\Omega)$ up to order $p$ (see [25]). For $X$ a Banach space, the Bochner spaces $L^{r}(0,T,X)$ are the spaces of functions $v:t\mapsto v(t)$ defined on $(0,T)$ with values in $X$ (see [25]). The functional space $\mathscr{W}_{r}$ and $\mathscr{M}$ are defined as

[TABLE]

with

[TABLE]

Regarding the mathematical analysis of this model, specifically the existence and regularity of solutions, we refer the readers to works such as [19, 20, 21] and the references therein. To insure existence of a solution of (1)-(3) in case of non–Newtonian fluids some assumptions on $\sigma({\boldsymbol{u}})$ are needed. The results proposed in [21], based on three general assumptions, is presented here. Since we are interested in generalized Newtonian fluid, two of these assumptions, here noted (H1) and (H2), are reformulated putting the emphasis on $\nu(t,{\boldsymbol{x}})$, the viscosity of the fluid. Notice that rheological models such as the power law and Carreau–Yasuda type models used in a wide variety of industrial applications all verify these assumptions with little restrictions (see Remark 1).

Hypothesis H1** (continuity)**

[TABLE]

is measurable with respect to $(t,{\boldsymbol{x}})$ for all $({\boldsymbol{u}},{\boldsymbol{F}})\in\mathbb{R}^{d}\times\mathbb{R}^{d\times d}$ and continuous with respect to $({\boldsymbol{u}},{\boldsymbol{F}})$ for almost every $(t,{\boldsymbol{x}})\in(0,T)\times\Omega$.

Hypothesis H2** (coercivity)**

There exist $c_{1}\geq 0,c_{2}>0,\lambda_{1}\in L^{r^{\prime}}((0,T)\times\Omega),\lambda_{2}\in L^{1}((0,T)\times\Omega),\lambda_{3}\in L^{(r/\alpha)^{\prime}}((0,T)\times\Omega),0<\alpha<r$, such that

[TABLE]

[TABLE]

with $\|{\bf F}\|^{2}=\sum_{ij}{\bf F}_{ij}^{2}$ the Frobenius norm of the tensor ${\bf F}$.

Hypothesis H3** (monotonicity)**

$\forall\ (t,{\boldsymbol{x}},{\boldsymbol{u}})\in(0,T)\times\Omega\times\mathbb{R}$* , the map ${\boldsymbol{F}}\mapsto\nu(t,{\boldsymbol{x}},{\boldsymbol{u}},{\boldsymbol{F}})$ is a $C^{1}$ function and is monotone,*

[TABLE]

Theorem 1

Assuming $\nu$ satisfies hypothesis (H1)-(H3) for some $r\in[\frac{3d+2}{d+2},\infty)$. Then for ${\boldsymbol{u}}_{0}\in\left(L^{2}(\Omega)\right)^{d}$ and ${\boldsymbol{f}}\in L^{r/(r-1)}(0,T;V_{r}^{\prime})$, the system (1)-(3) has a weak solution $({\boldsymbol{u}},p)$, with $u\in\mathscr{W}_{r}$ and $p\in\mathscr{M}$.

Remark 1

Following [21, 23, 22] the conditions (H1)-(H3) are verified in the case of a generic law

[TABLE]

Thanks to the numerous parameters, this correspond to various rheological models (power law, Carreau, Cross, Carreau–Yasuda, etc.) which are frequently used in engineering context. In such cases, from [22] we get existence of a weak solution, provided $m>0$ for $d=2$ and $m>1/5$ for $d=3$.

2.2 Time discretization

The choice of time discretization is motivated by the fact that we rely on splitting technique to construct approximation of (1). Projection schemes (as methods based on operator splitting) have an inherent splitting error of order $3/2$ in $H^{1}-$norm (see [26, 27, 1]). Therefore the proposed algorithm, relying on a projection scheme, is at best of second order in $H^{1}-$norm, and the use of higher order time discretization is irrelevant.

For the sake of simplicity and clarity, the method proposed in this paper will be based on the second order backward time discretization (BDF2) which is frequently used to solve Navier–Stokes equation (1).

[TABLE]

Using a constant time-step $\Delta t$, denoting $t^{n}=n\Delta t$, ${\boldsymbol{u}}^{n}={\boldsymbol{u}}(t^{n},{\boldsymbol{x}})$, $p^{n}=p(t^{n},{\boldsymbol{x}})$, ${\boldsymbol{f}}^{n}={\boldsymbol{f}}(t^{n},{\boldsymbol{x}})$ and $\nu^{n}({\boldsymbol{u}})=\nu(t^{n},{\boldsymbol{u}})$ the implicit time discretization of (1) result in a sequence of (generalized) Oseen problems of the form

[TABLE]

completed with the same initial and boundary conditions (2)-(3). Observe that since the viscosity could depend on the velocity (or pressure), (5) contains possibly two non linearities. As for the first time step $t^{1}$, a simple backward Euler time step could be used.

3 Toward a projection scheme for non–Newtonian fluid

In this section, we follow the approach used in [16] leading to the creation of a more precise projection scheme. This new scheme is a generalization of the incremental projection, as presented in [1] and its construction fellow the idea of the rotational projection in [2].

Starting with the incremental projection, at time $t^{n+1}$, the method consists in the following series of steps: solving the viscous nonlinear system (or Burgers’ equation) (6) gives a velocity prediction, next is a projection step (7) producing a pressure correction and finally the pressure updating step (8).

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

Remark 2

We could choose to replace $\nu^{n+1}({\boldsymbol{\tilde{{\boldsymbol{u}}}}})=\nu(t^{n+1},{\boldsymbol{\tilde{{\boldsymbol{u}}}}})$ by $\nu^{n+1}({\boldsymbol{u}}^{n})$ in (6), making the non linearity related to the viscosity disappear. Obviously using $\nu^{n+1}({\boldsymbol{u}}^{n+1})$ makes the approach almost unusable as it would impose to solve (6)–(7) as a coupled system.

The first objection regarding (6)–(8) is the fact that an artificial homogeneous Neumann boundary condition is enforced on the pressure. However, as illustrated in Figure 2–3, this inconsistent boundary condition as little effect on the overall precision of the scheme (see Rannacher et al. [8] for a detailed review of this point) and the $H^{1}$–norm of the error on the velocity prediction has a good behaviour.

More importantly this splitting neglect information easily at our disposal (related to the shear rate of the velocity prediction). The approximation resulting from (6)–(8) can be improved simply by recovering this information. In [2] the rotational pressure correction scheme is proposed. This scheme gives a quasi-consistent boundary condition for the pressure and better precision for the velocity and pressure, see [1]. However, this approach is only valid for homogeneous viscosity.

The more general shear rate projection proposed in [16], makes it possible to achieve this improved accuracy for heterogeneous viscosity. But as presented it does not take into account the explicit dependency of the viscosity upon velocity, leaving again valuable information that could enrich the pressure approximation. Assuming that the viscosity is heterogeneous and depends on the velocity field ${\boldsymbol{u}}$, we consider the new projection

[TABLE]

[TABLE]

[TABLE]

We will now establish the equation characterizing $\psi$ to obtain a splitting. From (10) and (11) we have

[TABLE]

Replacing in (9) we get

[TABLE]

then we want $\psi$ satisfying

[TABLE]

Taking the divergence, we have

[TABLE]

The right hand side of this equation is known, as for the boundary conditions for $\psi$, the compatibility condition for the Poisson equation gives us a Neumann condition (see (20)). For ${\boldsymbol{f}}$ sufficiently regular, using the Agmon-Douglis-Nirenberg theorem (see [28]), the solution of this Poisson problem, is well defined (see Remark  6) and we can establish that $\psi$ is of the same nature and as regular as the pressure.

Since the convective term is on ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}$, as in [2, 1, 16], from (9) and (12) we get a nearly consistent boundary condition for $p^{n+1}$ on $\Gamma_{D}$. A consistent condition would be obtained provided the non linear convective term in the momentum equation is neglected; which correspond to an unsteady Stokes problem.

[TABLE]

Remark 3

Obviously in case of solenoidal vector fields the right hand side of (13) could be further simplified. However, a priori, the velocity prediction ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}$ is not divergence free. From (13), $\psi$ is composed of two corrections related to the difference between ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}$ and ${\boldsymbol{u}}^{n+1}$: the first is associated with the variation of viscosity, the second one related to the variation of shear rate.

4 Strategies to solve the shear rate projection for non–Newtonian fluids

The prediction step (9) imply solving a problem containing two nonlinearities: the convective term and the viscous term (since $\nu$ depends on the velocity and possibly the pressure, etc.). Introducing ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}^{\star}$ and $\nu^{\star}$ to lighten the notation, (9) is written

[TABLE]

which encompass various strategies to solve (9) according to the definition given for ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}^{\star}$ and $\nu^{\star}$.

Defining ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}^{\star}$ and $\nu^{\star}$ as expression not depending on ${\boldsymbol{u}}^{n+1}$ would lead to ”explicit” strategies, conversely using expression depending on ${\boldsymbol{u}}^{n+1}$ or ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}$ would be called ”implicit” strategies. These implicit approaches correspond to some fixed point methods applied to solve (14): at time step $t^{n+1}$, given ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}_{0}$ an initial approximation of ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}$ solution of (9), until convergence, we seek ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}_{k+1}$ solution of

[TABLE]

If the choice of treatment can impact the overall order of convergence of the time discretization, it has no fundamental consequence on the viscous correction presented. We propose a few examples of treatment of both terms.

4.1 Treatment of convection term

Explicit approaches consists in using extrapolation (backward in time) for ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}^{\star}$, leading to a complete linearization of the convective term. The simplest approach consist in using, at time step $t^{n+1}$,

[TABLE]

A more effective and accurate strategy is the use of the Richardson extrapolation

[TABLE]

For implicit approaches an extrapolation (backward with respect to the iterative steps) of the convective term is used. The choice of ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}^{\star}$ will alter the behaviour of the iterative process. The simplest extrapolation consists, at time step $t^{n+1}$ for the iterative step $k+1$, in using

[TABLE]

Once again a Richardson extrapolation produce another strategy

[TABLE]

Finally, Newton’s method is treated separately as it would add terms to (15). In this case the use of ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}^{\star}$ can be dropped and (15) is replaced by (16): at time step $t^{n+1}$, given ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}_{0}$ an initial approximation of ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}$ solution of (9) and until convergence, we seek ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}_{k+1}$ solution of

[TABLE]

4.2 Treatment of the viscous term

The same basic strategies can be applied to the non linear viscous term. Explicit approaches consists in using extrapolation (backward in time) for $\nu^{\star}$. The simplest approach at time step $t^{n+1}$, is

[TABLE]

Two kind of Richardson extrapolations can be considered

[TABLE]

For implicit approaches, backward extrapolations (with respect to the iterative steps) of the viscosity are introduced. The choice of $\nu^{\star}$ will alter the behaviour of the iterative process. The simplest extrapolation at time step $t^{n+1}$ for the iterative step $k+1$, is

[TABLE]

Once again a Richardson extrapolation produce other strategies

[TABLE]

Finally, assuming we have sufficient regularity of $\nu$ with respect to ${\boldsymbol{u}}$, we can apply Newton’s method. Introducing

[TABLE]

the gradient of $\nu$ with respect to ${\boldsymbol{u}}$ at time $t^{n+1}$. Dropping the $\nu^{\star}$ notation and replacing (15) by (17) we get: given ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}_{0}$ an initial approximation of ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}$ solution of (9), we seek ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}_{k+1}$ solution of

[TABLE]

Both non linear terms can be treated explicitly, implicitly or by ”hybrid” strategy (using an explicit approach for one term and an implicit for the other). Notice that the use of explicit strategies for both terms makes (14) a linear equation, in all other cases we have a fixed point loop of general form (15) to consider (including a possible combined expression composed of (16) and (17)).

Remark 4

It is tempting to follow the same idea for the pressure, replacing $p^{n}$ in (9) and (11) by $p^{\star}=2p^{n}-p^{n-1}$ a second order extrapolation. However, as observed in [1], in that case the stability seems to necessitate a lower bound on $\Delta t$.

4.3 A shear rate projection algorithm for heterogeneous viscosity

We now summarize these different strategies in an algorithm, keeping in mind that certain choice for $\nu^{\star}$ and $\tilde{{\boldsymbol{u}}}^{\star}$ can modify the expression (18) in Step 2 below, or make it a linear system (in case of explicit expression). As a generic strategy for solving (9) we get the following algorithm: at each time step $t^{n+1}$,

Initialization:

[TABLE] 2. 2.

Solving the non linear equation in ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}^{n+1}$:

until convergence, compute ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}_{k+1}$ solution of:

[TABLE]

When converged put ${\boldsymbol{\tilde{{\boldsymbol{u}}}}}^{n+1}={\boldsymbol{\tilde{{\boldsymbol{u}}}}}_{k+1}$. 3. 3.

Projection step :

Compute $\varphi$ solution of Poisson problem:

[TABLE] 4. 4.

Velocity correction :

[TABLE] 5. 5.

Shear rate projection step :

Compute $\psi$ solution of Poisson problem:

[TABLE] 6. 6.

Pressure correction :

[TABLE]

Remark 5

For the first time step ($n=1$) a simple backward Euler or Crank-Nicolson scheme is used. Step 4 and 6 have to be interpreted as $L^{2}$-projections in the velocity and pressure space respectively.

Remark 6

For $\Gamma_{N}=\emptyset$ the Poisson problem in step 3 is ill-posed and we seek a solution in $H^{1}(\Omega)\backslash\mathbb{R}$. For $\Gamma_{N}\neq\emptyset$ the homogeneous Dirichlet boundary condition in (19) is imposed by the Helmoltz decomposition (as in [29]). From [4], step 5 has a unique solution in $H^{1}(\Omega)\backslash\mathbb{R}$ and its variational form is

[TABLE]

From a computational point of view, (19) and (20) can be viewed as solving the same problem twice with two different right hand side. For methods such as the finite element, strategies reducing the computational effort involved in solving both solutions are available.

Remark 7

The basic incremental projection [14], consists in ignoring step 5 and putting $\psi\equiv 0$ which gives the original system (6)-(7). ²For an homogeneous viscosity, solving (20) is not required as in this case, thanks to the properties of the rotational operator, we have

[TABLE]

*which is the rotational projection [2]. In case of a heterogeneous viscosity, not depending on ${\boldsymbol{u}}$, obviously the first term in (13) will cancel, leaving only the second term resulting in the projection presented in [16]. Finally, in step 6 (pressure correction), applying a ”relaxation” factor $\alpha\in]0,1]$ to $\psi$ we get a generalized version of the ”Chorin-Uzawa” scheme as described by Rannacher in [26]. *

5 Finite element discretization

For the spatial approximation, the finite element method is used. First, we introduce a mesh $\mathcal{T}_{h}=\{\top\}$ of simplicial $\top$ partitioning $\Omega$:

[TABLE]

with the usual restrictions on $\mathcal{T}_{h}$ (see [30] or [31] for details). We emphasize that at no point in what follows are we making assumptions on the nature of the mesh. This partitioning of $\Omega$ could be composed of triangles, quadrangles or both for $\Omega\in\mathbb{R}^{2}$ or tetrahedra, hexahedra or both for $\Omega\in\mathbb{R}^{3}$. From the family $\mathcal{T}_{h}$ of partition of the domain $\Omega$ (indexed by $h$), we construct the family of finite dimensional vector spaces.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $P_{k_{\alpha}}(\top)$ is the space of polynomials of degree less or equal to $k_{\alpha}$ on $\top$. The index ”$h$” is used to denote the different spatial approximations. The approximation of $f^{n}({\boldsymbol{x}})$, a function $f(t,{\boldsymbol{x}})$ at time $t^{n}$, will be denoted $\text{f}^{n}_{h}$. At each time step $t^{n+1}$, we are now seeking approximations

[TABLE]

Since the discrete inf–sup condition must be respected for $({\bf u}_{h},\text{p}_{h})$ [4, 1, 3], the prediction equation (18) (recall that (16), (17) or a composite expression from both could replace (18)) is discretized using specific choice of interpolations, typically $k_{p}=k_{\bf u}-1$, see [3] for the details and other choices. Here a second order Taylor-Hood $P_{2}-P_{1}$ element was chosen for ${\bf V}_{h}$ and $M_{h}$. Regarding the approximation of $\varphi$ and $\psi$, a priori, they could be of arbitrary interpolating degree. Even though steps 4 and 6 are interpreted as projections in ${\bf V}_{h}$ and $M_{h}$, we must take into account that $\varphi$ and $\psi$ are solutions of Poisson problems and that they act as velocity and pressure corrections. In this work, $\varphi_{h}$ and $\psi_{h}$ were chosen to be quadratic ($k_{\varphi}=k_{\psi}=2$).

Following the usual approach, the Galerkin formulation based on ${\bf V}_{h}$, $M_{h}$, $M_{h}^{\varphi}$ and $M_{h}^{\psi}$ is used to generate finite dimensional systems corresponding to (18), (19) and (20). As underlined in Remark 6, the discrete version of Step 3 and Step 5 must be treated in a special manner (at least one discrete system is under-determined) and a unique solution is obtained through a supplementary computation in the last step (see Step 6h). To complete the spatial discretization of the shear rate algorithm a proper interpretation of Step 4 and Step 6 must be given. Here we chose the most natural approach which consist in the use of a $L^{2}$–projection of $\varphi_{h}$ and $\psi_{h}$ in $M_{h}$ (Step 6) and of $\nabla\varphi_{h}$ in ${\bf V}_{h}$ (Step 4) (see Remark 9). We are now in measure of giving a completely discrete generic version of the shear rate projection algorithm:

Initialization:

[TABLE] 2. 2.

Solving the non linear prediction equation:

until convergence, compute ${\tilde{\bf u}}_{k+1}\in{\bf V}_{h}$ solution of

[TABLE]

When converged put ${\tilde{\bf u}}_{h}^{n+1}={\tilde{\bf u}}_{k+1}$. 3. 3.

Projection step :

Compute $\varphi_{h}\in M_{h}^{\varphi}$ solution of

[TABLE] 4. 4.

Velocity correction :

Compute ${\bf u}^{n+1}_{h}\in{\bf V}_{h}$ solution of :

[TABLE] 5. 5.

Shear rate projection step :

Compute $\psi_{h}\in M_{h}^{\psi}$ solution of

[TABLE]

with ${\bf S}_{h}=\left(2\nu^{\star}{\boldsymbol{D}}({\tilde{\bf u}}_{h}^{n+1})-2\nu^{n+1}({\bf u}^{n+1}_{h}){\boldsymbol{D}}({\bf u}^{n+1}_{h})\right)$ in $\left(M_{h}\right)^{d\times d}.$ 6. 6.

Pressure correction :

Compute $\text{p}^{n+1}_{h}\in M_{h}$ solution of

[TABLE]

with

[TABLE]

Remark 8

For $n=1$, the use of a backward Euler or Crank-Nicolson scheme slightly modifies the expressions in Step 2h, 3h and 4h. In Step 5h, ${\bf S}_{h}$ is defined through an $L^{2}$ projection. In Step 6h, $\mu(\cdot,\cdot)$ correspond to an average on $\Omega_{h}$ of $\psi_{h}$ or $\varphi_{h}+\psi_{h}$. This correction gives $\psi_{h}\in M_{h}^{\psi}\backslash\mathbb{R}$ (and $\varphi_{h}\in M_{h}^{\varphi}\backslash\mathbb{R}$ if needed) insuring the uniqueness of the projection in Step 6.

Remark 9

Concerning the sum in Step 6h, a possible alternative, avoiding solving an algebraic system, would be to employ a hierarchical basis for $M_{h}^{\varphi}$ and $M_{h}^{\psi}$ and to use the $P^{1}$ component of $\varphi_{h}$ and $\psi_{h}$ instead of a $L^{2}$–projection. Such approach, if tempting considering computational costs, could result in important loss of information, negating the effort made in determining $\psi_{h}$ at Step 5.

6 Numerical experiments

The goal of this section is to validate this family of methods numerically, illustrate the accuracy of some representatives and compare them. Two tests will be presented:

the first is an analytic test, based on a two dimensional manufactured solution presenting the time and spatial accuracy of four representive methods.

• 2.

The second test is the simulation of a flow past a cylinder of a Carreau fluid for various values of the power index [32]. Here the drag and wake (or re-circulation) lenght will be the object of comparison.

6.1 Choice of representatives and efficiency

Section 4 offers various approaches for the shear rate projection. These different strategies influences accuracy, efficiency and computational costs. It is not our intention to establish the relative merit of all those methods. We will merely illustrate the capability of a small number of schemes when used on specific problems. For the first test we chose fully explicit and fully implicit representatives, leaving hybrid and Newton’s variation to future, more thorough, analysis. For the second test, as we try to reproduce a steady state, the unsteady implicit shear rate projection is used and compared to a ”steady” mixed algorithm and to results from the literature.

The use of the explicit strategy for both non linear terms makes (14) and (21) linear equations. This is tempting as it reduces computational efforts at each time step. However, totally explicit (as well as semi–explicit (hybrid) strategies) leads to conditional stability and certainly imposes limitations on the lenght of the time step. Nevertheless it would be erroneous to exclude such approaches, as they can translate into very efficient algorithms.

Implicit strategies, fixed-point or Newton-like approaches, are of general use. Their efficiency is sometimes questioned since these choices leads to an iterative process. However implicit discretization minimizes the constraint on time step length making it possible to regain (when compared to explicit methods) numerical efficiency through the use of larger timestep (for a comparable precision). In those cases, it results in a better overall performance of the algorithm.

Four schemes will be used,

a basic explicit incremental projection, denoted (IPex), obtained by imposing $\psi_{h}\equiv 0$, ignoring Step 5h and putting, at time $t^{n+1}$

[TABLE] 2. 2.

a basic implicit incremental projection, denoted (IPim), obtained by imposing $\psi_{h}\equiv 0$, ignoring Step 5h and putting, at time $t^{n+1}$

[TABLE] 3. 3.

the explicit shear rate projection scheme, denoted (SRPex), obtained by imposing, at time $t^{n+1}$

[TABLE] 4. 4.

the implicit shear rate projection scheme, denoted (SRPim), obtained by imposing, at iteration $k+1$ of time $t^{n+1}$

[TABLE]

The algorithms are implemented in a finite element context using the FreeFem++ software, see [24]. Note that all computations are performed using triangular mesh and Taylor-Hood interpolation: velocity are quadratic interpolation ($P_{2}$) and the pressure is linear ($P_{1}$). Insuring the respect of the inf–sup condition (see [1, 3]). Both projections, $\varphi$ and $\psi$, are taken in $P_{2}$.

6.2 Accuracy test

The first test is used to illustrate and compare the time and space accuracy of the shear rate projection (SRP) with the basic incremental projection (IP) for viscosity depending on the velocity of the fluid. As in [16] a variation on the finite element tests proposed by Guermond et al. in [27] is proposed

[TABLE]

with a suitable forcing term depending on $\nu$.

6.2.1 Time accuracy test

We use a uniform mesh sufficiently fine, $200\times 200$, to insure negligible spatial error for the chosen range of time step. As for Step 2h, when an implicit form is used, the stopping criterion of the fixed point uses the norm of the variations of ${\tilde{\bf u}}$ with a tolerance of $10^{-8}$. The velocity and pressure error to the exact solution will be computed using two types of norms, denoting $e$ the error, we will use

[TABLE]

As expected, see Figure 2–3, the explicit form of both methods, IPex and SRPex, exhibit a loss of accuracy when compared to their implicit version. The pressure error, Figure 1 and 3 illustrate the gain in accuracy of this ”extended” SRP method, this is comparable to the results presented in [16].

For the velocity in $L^{2}$–norm, (Figure 2 right), both explicit methods exhibit a linear rate while both implicit methods have a quadratic convergence rate. The SRP does not produce any significant gain in precision for the velocity. The $H^{1}$ convergence rates are identical: the explicit methods have a rate of $1.1$ while the implicit method shows a rate of almost $2$, this is explained by the increased accuracy of the pressure. The SRPim produces a velocity error in $H^{1}$ norm approximately $2$ time smaller then the $H^{1}$ error produced by the implicit incremental projection.

The $L^{2}$ rates for the pressure (Figure 3 left) seems relatively similar for both explicit and implicit pair of schemes (around 1 for the explicit schemes and 1.7 for the implicit ones). However both version of the SRP are clearly more accurate, reducing the $L^{2}$–error by a factor roughly equal to $3$. From Figure 3 the $L^{\infty}$ rate of both IP methods is linear and the SRP’s rate are $3/2$: the extreme values of the errors of the SRP are decaying more rapidly then those of the IP method. Figure 1 illustrate at time $t=1$, for $\Delta t=1/80$, the spatial distribution of the error for the pressure; giving a sense of the impact of $\psi$ on the pressure.

In conclusion, for non-homogeneous viscosity, both version of the SRP method produce a more accurate approximation of the pressure (and velocity) than the IP version, even in the explicit case, the use of the SRP improves the pressure approximation.

6.2.2 Spatial accuracy test

Since the SRP schemes contain an additional problem to solve, (23), it is necessary to illustrate the impact of this added computations on the spatial accuracy of the scheme. Here we included two other methods as comparative: IPex which is the simplest transition to projection methods for generalized Newtonian fluids and which has the smallest computational cost per time step of all the methods presented and an implicit version of the very simple penalty method [7], also called mixed method, using a penalty of $10^{-10}$ for the incompressibility.

For this accuracy test, 6 meshes obtained by regular subdivision from $1/10$ to $1/320$ were used. A total of 50 time steps of fixed lenght $10^{-4}$ (for a final time t=0.005) was used; assuming such small time step would limit the impact of the time discretization on the approximation error. The first observation concerns the IPex scheme, the $L^{2}$–error for the velocity and pressure seems to be rapidly saturated, mesh refinement over $1/80$ having no effect on the error for the current time step of $10^{-4}$. From the finest mesh, the SRPim seems to exhibit a faster saturation of the $L^{2}$–error compared to the penalty method (denoted Mixed in Figure 4). Nevertheless, the SRPim and Mixed method produce solutions with comparable spatial convergence rate, confirming the good behaviour in space and in time of $\psi_{h}$ solution of (23).

6.3 Flow past a cylinder

For this test, a two-dimensional model of the flow of a generalized Newtonian (Carreau) fluid past a cylinder, see figure 5, is considered. This second test uses more realistic data (boundary conditions, viscosity, etc.) and can be viewed as a ”reference test” as it is presented frequently (see [32, 33, 34] and the references therein, or [35] for a more complex law (Herschel-Bulkley)).

This well documented test offers reference values for punctual and boundary related physical measures. This allows us to illustrate the accuracy of the projections schemes based on relevant ”local” physical quantities of interest: the drag coefficient $C_{D}$ and the wake lenght $L$,

[TABLE]

The fluid is characterized by the Carreau law, defining the viscosity, and by two adimensional numbers, the Renolds number ($Re$) and the Carreau number $C_{U}$,

[TABLE]

For this test we are interested in a steady state flow, which correspond to a relatively small $Re$. Contrary to tests frequently encountered in the literature, we will limit our tests to a single Reynolds number. Experimentation will be made using the implicit shear rate projection (showing the ability of the SRPim to reproduce steady states results). Reference solutions will be computed using a steady state mixed algorithm (denoted SS-Mix) based on a penalty method inside a fixed point loop for the viscosity.

The data for the tests are taken from [32, 34], adding another point of comparison. We can summarize the parameters as follow:

[TABLE]

Finally, contrary to  [32, 33], a natural (or open) boundary condition is applied on $\Gamma_{o}$.

As an unconfined flow is modeled, and considering the open boundary condition, the outer boundary of the domain need to be sufficiently far from the cylinder so that boundary effects are at a minimum. This implies a relatively large domain. It should be notice that the dimensions chosen guarantees us of eliminating any boundary effects; however the use of a smaller domain could be possible (provided a proper domain dependency analysis is made). Mesh adaptation is used to reduce the number of elements by refining and de-refining locally the mesh, this makes it possible to have reasonable computational time even for such large domain. For each power index value ($m$), a specifically adapted mesh will be used. These meshes are constructed based on the SS-Mix solutions for each value of $m$ (see [24] for details concerning the mesh adaptation procedure).

From [32], the critical Reynolds, at which the flow becomes unsteady, for this Carreau fluid and for a power index $m\geq 0.4$ is greater then $10$. Therefore the flows simulated here are certainly steady for all the values of $m$ used, the steady state will be reached using the SRPim scheme, an unsteady scheme. The simulation is stopped once the steady state is attained, as indicated by the stationnarity of certain physical quantities (here the drag coefficient and magnitude of the velocity).

Figure 6 presentes some streamlines, near the obstacle, for the velocity obtained with the SRPim for both $C_{U}$ number, with $m=0.4$ and $0.6$. This illustrates the variations of the wake lenght and overall behaviour of the method which is in agreement with the literature of the subject.

Table 1 presents the drag coefficient and wake (or re-circulation) lenght for both values of the Carreau number ($C_{U}=10$ and $20$) and a power index ranging from $0.4$ to $1.0$. Relative variations of the results from the steady state mixed method and the SRPim are also presented. In both cases these variations are roughly less than $1\%$, which can be attributed to the convergence criteria to reach the steady state, confirming the good behaviour of the SRPim method even for local quantities of interests such as the drag coefficient. Although not presented here, the results of Table 1 for $C_{U}=10$ and $m=0.4,0.6,0.8,1.0$ are in good agreement with the results presented in [32] and [34] with less then $2.5\%$ of variation for the drag coefficient.

7 Conclusions

We proposed an original projection method for the numerical simulation of generalized Newtonian flow. This method takes into account the explicit dependence of the viscosity upon velocity. Based on this projection and the finite element method, a family of numerical scheme using a second order time approximation was proposed.

Using a manufactured solution and the simulation of the flow past a cylinder, we illustrate the validity of four representatives of the family of numerical methods. As expected, explicit methods were less expensive from a computational point of view, however they exhibit the usual weakness of such methods (mainly weak order of precision). As for the implicit versions, the SRPim can be viewed as a generalized rotational projections. The behaviour of both IPim and SRPim was satisfactory. We underline that the SRP exhibit, as expected from [16], an improved accuracy in approximating the pressure when compared to the more intuitive incremental approach.

Finally, as in [36], these shear rate projections could be used to produce new coupled projection methods applicable to generic convection-diffusion phenomenon coupled with Newtonian/Non Newtonian flows (natural convection, coextrusion, etc.). Following the analysis in [36] these methods should offer important gain in performance and accuracy.