SciELO - Scientific Electronic Library Online

vol.22 número1Generalized line criterion for Gauss-Seidel methodApproximate controllability for the semilinear heat equation in R N involving gradient terms índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Computational & Applied Mathematics

versión On-line ISSN 1807-0302

Comput. Appl. Math. v.22 n.1 São Carlos  2003


Finite element approximation of bipolar viscous fluids



H. ManouziI; A. BrahmiII; M. FarhloulII

IDépartement de Mathématiques et de Statistique Université Laval, Québec, G1K7P4, Canada
IIDépartement de Mathématiques et de Statistique Université de Moncton, Moncton, NB, E1A 3E9, Canada




A bipolar viscous fluid model is assumed to regularise the solution of Newtonian and quasi-Newtonian flows. In this article, a mixed finite element approximation of the bipolar viscous fluids is proposed. In this approximation the velocity of the fluid together with its laplacian are the most relevant unknowns. An existence and uniqueness results are proved. A mixed finite element approximation is derived and numerical results are presented.

Mathematical subject classification: 65N30, 76M10.

Key words: finite element, bipolar fluids, Navier-Stokes, mixed finite element, PDE.



1 Introduction

Mathematical model for fluid motion play an important role in theoretical and computational studies in the aeronautics, meteorological, plasma physics, etc. The Navier-Stokes equations are generally accepted as providing an accurate model for the incompressible motion of viscous fluids in many practical situations. In the classical Navier-Stokes equations theory, the viscosity of the fluid is modeled by the dependence of the stress on the first gradient of velocity. It is well known that the corresponding mathematical theory contains a number of unsolved problems. One of the fundamental open problems for the Navier-Stokes equations in 3-dimensional space is that of the uniqueness of the solution of the Cauchy-Dirichlet problem, which is not guaranteed in the functional classes in which global existence holds. These problems have led many authors that a stronger mechanism of dissipation and viscosity, namely the dependence of the stress on the higher gradients of velocity, must occur in the flows of viscous fluids. In materials in which the higher gradients of velocity and the higher gradients of deformation influence the response, the rate of work of internal forces cannot be expected to be only the product of the usual second order stress tensor with the first gradient of velocity. Instead of this, a more general expression must be assumed containing additionally the sum of products of higher order multipolar stress tensors with the higher gradients of velocity. Otherwise such materials cannot be compatible with the Clausius-Duhem inequality. In [13] a thermodynamics theory of constitutive equations of multipolar viscous fluids within the framework of the theory of Green and Rivlin [9] has been developed. This paper is organized as follows. In section 2 we formulate the problem, we introduce some notations and we give some preliminary results. In section 3 we formulate the bipolar Navier-Stokes problem and we prove the existence and the unicity of the solution. In section 4 we study a mixed variational formulation and its approximation. Finally in section 5 we give some numerical results.


2 Problem formulation and preliminary results

2.1 Problem formulation

We will use the standard notation for gradient, divergence and curl operators, i.e. for a scalar function p and a vector function v, we write Ñq, div v and curl v for the gradient, divergence and the curl respectively. We further recall that the vector-valued Laplace operator, D, is defined by

Dv = (Dv1, Dv2), v = (v1, v2)

and we have the vector identity

Next we let W be an open bounded domain of 2, with smooth boundary G. We consider the stationary flow of an incompressible nonlinear bipolar fluid. This problem is modeled by the following system of equations:

where u is the velocity field, s the stress tensor and f represents the external forces. The stress tensor s = (sij) has the form (see [4], [8], [9],[10], [11], [13])


and r, n, l and m are positive constants with r > 1. Since div u = 0, the divergence of the stress tensor takes the form

divs = -Ñp + 2n div (F r((u)) (u) ) - mD2u.

Correspondingly, the bipolar Navier-Stokes equations become

By choosing the appropriate constants r, m, n and l we find these classical models of fluids:

  • Carreau's law model (see [15]): m = 0.
  • Power law model (see [15]): m = 0, l = 0.
  • Classical Navier-Stokes: r = 2, m = 0.
  • Regularized Navier-Stokes equations (see [14]): r = 2 , m ¹ 0.

These laws occur in many mathematical models of physical processes of polymer, of visco-elastic and visco-plastic fluids or of fluids with large stresses (See [10]). We will study the system of equations (4) with the boundary conditions

where t is the unit tangent to G.

2.2 Preliminary results

Let Hm(W) and Hm(G) denote the usual Sobolev spaces of order m defined on W and G, respectively; also Hm(W) and Hm(G) denote their vector-valued countreparts. The inner product on Hm(W) or Hm(W) is defined by (., .)m. The space Hm(W) is equipped with the norm


For vector valued functions the above norm extends naturally as follows: for v = (v1, v2)

We also define the Hilbert spaces

H (div, W) = {v | v Î L2(W), div v Î L2(W)}
H (curl, W) = {v | v Î L2(W), curl v Î L2(W)}

These spaces are equipped with their natural hilbertian norms.

Let (W) denote the space of infinitely differentiable functions having compact support and let H(W) be the closure of (W) with respect to the norm . m,W and Hm(W) the dual space of H(W). The duality pairing between H(G), H(G) and their dual spaces H(G), H(G) is denoted by á., .ñ , , G.

We also have the following trace results (see, e. g., [6], [16])

Theorem 2.1. Let W be a bounded open set of 2 with a Lipschitz continuous boundary G. We have the following properties:

1. There exists a unique linear continuous map g0 : H1(W) H(G) such that

g0v = v|G, "v Î H1(W).

2. There exists a unique linear continuous map gn : H(div, W) H(G) such that

gn v = (v · n)|G, "v Î H(div, W)

where n denotes the unit outward normal to G.

3. There exists a unique linear continuous map gt: H(curl, W) H(G) such that

gt v = (v · t)|G, "v Î H(curl, W).

We use the convention H0(W) = L2(W), H0(W) = L2(W), (., .)0 = (., .).

Also we introduce some function spaces which are related to the study of our problem:

For details concerning Sobolev spaces consult [1], [5], [6].

Lemma 2.1. Let u Î L2(W) such that Du Î , then we have the following Green Formula:

Proof. The proof of this equality follows from the classical Green formula:

and from the fact that curl curl v = –Dv and curl curl Du = –D2u for all u and v satisfying the conditions of lemma 2.1.

The following preliminary results are classical and will be assumed ( see, e.g., [5], [7]).

Lemma 2.2. Let W be bounded with a smooth boundary G. Then the mapping

v Dv0,W

defines a norm on X which is equivalent to the norm induced by H2(W), i.e.

$a > 0, Dv0,W > avX.

From lemma 2.2 we see that in X the norm u2,W induced by H2(W) is equivalent to the norm Du0,W. We then denote

as a norm in X derived from the inner product (u,v)X = (Du, Dv)0,W.

Lemma 2.3. Let W be bounded with a smooth boundary G. Then there exists a positive constant C = C(W) such that


3 The bipolar Navier-Stokes problem

3.1 Variational formulation

A variational formulation of the boundary value problem (4) and (5) reads as follows:

Remark 3.1. The term

is well defined for every (u,v) Î H2(W) × H2(W).

We associate with the problem (7) the following problem

Theorem 3.1. We have the following results:

1. If (u, p) is a smooth solution of (4)-(5) then u is a solution of (8).

2. If u is a smooth solution of (8) then there exists a function p Î M such that (u, p) is a solution of (4)-(5).

Proof. 1. Let v Î V, if (u, p) Î X × M is a solution of (4)-(5) with curl Du Î H1(W) then we have

We derive from the Green formula (lemma 2.1) that

Since v Î V, div u = 0 and Du · t = 0, we have u Î V and

This proves the first part of the Theorem.

2. Now, let u Î V be a solution of (8) with curl Du Î H1(W) then, we have

Let us introduce the operator

Note that < L, v > = 0 "v Î and, since -(Du, Dv) = (Du, curl curl v) = (curl Du, curl v), we get L Î H1(W). Then, since W is simply-connected and from De Rham's lemma, there exists one and only one function p Î M such that

L = –Ñ p.

Then, we readly derive that

Now, Owing to Green's formula (lemma 2.1), we have

Hence, it follows that

Since f Î L2(W), we have

Moreover, as u Î V, we have

Now, by using the Green's formula (lemma 2.1), we have on the one hand:

and, by (8), we get

Finally, from [2], we have

gt(Du) = 0 on G.

3.2 Existence and unicity of the solution

In order to show the existence and the unicity of the the solution of problem (8) we need to state some intermediate results.

Definition 3.1. Let X be a reflexive Banach space and let A be an operator from X to X¢. We say that

  • The operator A is monotone if and only if

< Au Av, u v > > 0 "(u, v) Î X × X

  • The operator A is pseudo-monotone if A satisfies (un u) and

(limp sup áAun, un – uñ < 0)
áAu, u – wñ < limp sup áAun, un – wñ

for all w Î X.

  • The operator A is strongly continuous if for every sequence unwith un u we have Aun Au.

Lemma 3.1. Let A1 and A2 be two operators from X to X¢ with

  • A1 is monotone
  • A2 is strongly continuous

Then, A1 + A2 is pseudo-monotone.

Proof. See [12], [17].

Lemma 3.2. Let X be a reflexive and separable Banach space and let A be an operator from X to X¢ with

  • A is continuous
  • A is coercive i.e.
  • A is pseudo-monotone.

Then, for every b Î X¢ the equation Au = b has at least one solution u Î X.

Proof. See [12], [17].

In the following lemma the notation (s, t) means .

Lemma 3.3. For any s, t Î 2×2, there exist constants C1, C2, C3, C4 such that:

1. If 1 < r < 2, then

2. If r > 2, then

Proof. See [15].

Lemma 3.4. Let u, v Î H1(W) and let w Î H1(W) with div w = 0 and (w · n)|G = 0. Then, we have

  • ((w · Ñ) \u, v) + ((w · Ñ)v, u) = 0
  • ((w · Ñ)v, v) = 0

Proof. See [6], [16].

Now, we introduce the operators A0, A1 and A2 defined, for all u, v Î V, by

  • (A0u, v) = m(Du, Dv)
  • (A1u, v) = ((u · Ñ)u, v)
  • (A2u, v) = 2n(Fr(e(u))e(u), e(v)).

Let the operator A = A0 + A1 + A2. Then, the problem (8) can be written as follows

Proposition 3.1. The operator A is pseudo-monotone.

Proof. According to lemma 3.1, we first prove that A1 is strongly continuous. To this end, let un be a sequence of elements of V such that un u in V. Owing to lemma 3.4, we have, for all v Î V,

(A1un, v) – (A1u, v) = (((un u) · Ñ)un, v) + ((u · Ñ)v, u – un)

and then

|(A1un, v) – (A1u, v)|
< C(un u1,Wun1,Wv1,W + u1,Wv1,Wu un1,W).

Now, since the embedding of V in H1(W) is compact, we have un u in H1(W). Therefore

(A1un, v) (A1u, v), "v Î V.

On the other hand, owing to the lemma 3.3, we prove that A0 + A2 is monotone. Indeed, for all u, v Î X we have

(A2u – A2v, uv) = (Fr(e(u))e(u) – Fr(e(v))e(v), e(u) – e(v)) > 0


(A0u – A0v, uv) = m(D(uv), D(uv)) > 0.

Therefore, owing to lemma 3.1, the operator (A0 + A2) + A1 i.e. A is pseudo-monotone.

Proposition 3.2. The operator A is continuous and coercive.

Proof. It is obvious, since A0, A1 and A2 are continuous, that A is continuous. Beside that, since ((v · Ñ)v, v) = 0, "v Î V, we have, for all v Î V,

In virtue of lemma 3.3, we have (Fr(e(v))e(v), e(v)) > 0. Then, we derive that

Finally, by using lemma 2.2 we conclude that A is coercive.

Theorem 3.2. The problem (8) has at least one solution u.

Proof. From propositions 3.1 and 3.2 we have

  • A is continuous
  • A is pseudo-monotone
  • A is coercive.

Then, in virtue of lemma 3.2 and since V is reflexive and separable, we have the existence, of u Î V such that Au = f, or equivalently the problem(8) has at least one solution.

Theorem 3.3. If m is sufficiently large or f ''sufficiently small" so that

m2 > C(W)fX¢

then problem (8) has exactly one and only one solution.

Proof. We assume that there are two solutions u1 and u2 in V of (8) then we have

(Au1 – Au2, v) = 0 "v Î V.

Then, for all v Î V, we have

m(Du1Du2, Dv) + 2n(Fr(e(u1))e(u1) – (Fr(e(u2))e(u2),
e(v)) + ((u1 · Ñ)u1, v) – ((u2 · Ñ)u2, v) = 0.

By taking v = u1u2, we get

Since, owing to lemma 3.3,

2n(Fr(e(u1))e(u1) – Fr(e(u2))e(u2), e(u1) – e(u2)) > 0,

we readily derive that

Now, for some positive constant C1 depending on W, we have

By using lemma 2.2, we get

On the other hand, since u1 is a solution of (8) and by taking v = u1 in (8), we get


We easily derive that u1 satisfies

and then

As a consequence of (25) and (26), we obtain


Finally, we readily derive the uniqueness of the solution when the condition



4 The mixed finite element formulation

In order to approximate numerically problem (7) we introduce the mixed formulation which allows us to reduce the order of problem (7). That makes the numerical approximation easier. Let us take w = –Du as a new variable. Then, since div (Du) = 0 and owing to the identity (1), we may recast problem (4) as second order system for u and w:

In order to study the last problem we introduce the following space

Q = H0(curl, W) = {v Î H(curl, W) and (v · t)|G = 0}.

A variational formulation of the mixed problem reads as follows

Theorem 4.1. 1) If (u, p) Î X × M is a solution of (4)-(5) and w = –Du with curl w Î H1(W) then (u, p, w) is a solution of problem (28).

2) If (u, p, w) is a solution of problem (28) then (u, p) is a solution of (4)-(5).

Proof. 1. Let (u, p) a solution of (4)-(5) and w = –Du with curl w Î H1(W) . Then, as in the proof of theorem 3.1, we have Du · t = 0 i.e. w · t = 0. This yields w Î Q.

Beside that, for all v Î , we have


Now, by substitution in (29) we get

Owing to De Rham's lemma there exists a unique function p Î M such that

On the other hand by using the new variable w = –Du and the fact that div u = 0, we get for all z Î Q

Then, taking into account that z · t = 0 and u Î X, we find

(w, z)-(curl u, curl z) = 0 "z Î Q.

Finally, from div u = 0, we get

(div u, q) = 0 "q Î M.

2. Conversely, suppose that (u, p, w) is a solution of (28) such that u Î H2(W) and curl Du Î H1(W). First, since (div u, q) = 0, "q Î M and u Î H(W), we have div u = 0. Now, since div u = 0, we have

Thus w = –Du, and Du Î Q since w Î Q. Therefore,

Thus, by the fact that div Du = Ddiv u = 0 , Du Î Q, curl Du Î H1(W) and by using the lemma 2.1, we have

(curl w, curl v) = (D2u, v) "v Î .


Owing to De Rham's lemma and like in the proof of theorem 3.1, there exists one and only one function p Î M such that:

Finally, by gathering, we get

Problem (28) is easier to approximate numerically than problem (7) since the finite element approximation of (28) involves only the construction of finite-dimensional sub-spaces of the spaces H0(curl, W) , H(W) and L2(W). In this section we will only consider the discrete approximations to the mixed variational formulation (28). We assume for simplicity that the domain W is a two-dimensional polygon. Let h be a given triangulation of into closed triangles, with h = max{diam (K)} where the maximumis taken over all K Î h. Let Pk(D) denote the space of polynomials of degree less than or equal to k with respect to the set D Ì 2 and Pk(D) denotes the space of 2-vector valued functions each of whose components belongs to Pk(D). To approximate the velocity and the pressure, we use the Crouzeix-Raviart finite element:

where, for any K Î h,

Here, l i(x), i = 1,2,3, denote the barycentric coordinates of the point x Î K with respect to the vertices of K. Finally we introduce the laplacian's velocity finite element space:

Now, let

Then (28) is discretized by the following problem


5 Numerical results

The performance of the numerical method described above has been tested on three examples of a comparative nature: the Poiseuille flow, the flow over a backward-facing step and the driven cavity flow.

  • Poiseuille flow

    At inlet of the channel, a parabolic profile for the velocity is used. Figures 1, 2 give the velocity fields for the Reynolds number Re = 1000 and for r = 1.7 and r = 1 respectively. The velocity distribution in the middle of the channel is shown in Figure 3 for different values of the parameter r = 1.8, 1.6, 1.4, 1.2, 1.1, 1.







  • Backward-facing step flow

    With flow in channels of constant geometry such as the Poiseuille flow, a predominant flow direction prevails. Instead, where the geometry changes abruptly, the flow separates and develops a recirculation region. Inlet value for the velocity is the same as that of Poiseuille flow. The Reynolds number is based on step height and a ratio outlet height to step height of 2. The streamline patterns are presented in Figures 4-11. The computations were performed for differents Reynolds numbers and for differents parameters r and m.

















  • The driven cavity flow

    The depicted domain is the square cavity [0,1] × [0,1], and the boundary conditions are, u · t = 1 and u · n = 0 at y = 1, and u = 0 on the other parts of the boundary. In Figures 12-17 the streamlines contours are given for different Reynolds number and different parameters r and m.













For all the above numerical tests the value of the parameter l has been fixed to the value l = 10-11.



[1] R.A. Adams, Sobolev spaces, Academic Press, New York, (1975).         [ Links ]

[2] C. Amrouche and V. Girault, Problèmes généralisés de Stokes, Publication du Laboratoire d'Analyse Numérique R91001, Université Pierre et Marie Curie, Paris.         [ Links ]

[3] H. Bellout, R. Bloom and J. Necas, Phenomenological behavior of multipolar viscous fluids, Quarterly of Applied Mathematics, Volume IX, Number 4 (1991), 687-728.         [ Links ]

[4] A. Brahmi, Simulation numérique des fluides bipolaires, Mémoire de maîtrise, Université Laval, Québec, (1994).         [ Links ]

[5] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978.         [ Links ]

[6] V. Girault and P.A. Raviart, Finite element for Navier-Stokes equations, theory and algorihms, Springer Ser. Comput. Math., 5, Springer-Verlag, Berlin, New Tork, (1986).         [ Links ]

[7] R. Glowinski and O. Pironneau, Numerical methods for the first biharmonic equation and for the two-dimensional Stokes problem, SIAM Rev., 21, 2 (1979), 167-212.         [ Links ]

[8] K.K. Golovkin, New Equations Modeling the Motion of a Viscous Fluid, and their Unique Solvability, Proc. Steklov Inst. Math., 102 (1967).         [ Links ]

[9] A.E. Green and R.S. Rivlin, Simple force and stress multipoles, Arch. Rat. Mech. Anal., 17 (1964), 325-353.         [ Links ]

[10] O.A. Ladyzhenskaya, Modification of the Navier Stokes Equations for Large Velocity Gradients, in Boundary Value Problems of Mathematical Physics and Related Aspects of Function Theory II, Consultants Bureau, New York, (1970).         [ Links ]

[11] O.A. Ladyzhenskaya, Nonstationary Navier-Stokes equations, Amer. Math. Soc. Transl. (2) 25 (1962), 151-160.         [ Links ]

[12] J.L. Lions, Quelques méthodes de résolution des problèmes aux limites non linéaires, Dunod, Paris, (1969).         [ Links ]

[13] J. Necas and M. Silhavy, Multipolar viscous fluids, Quart. Appl. Math., 49 (1991), 247-265.         [ Links ]

[14] Y.U. Ou and S.S. Sritharan, Analysis of regularized Navier-Stokes equations I, Quart. Appl. Math., 49 (1991), 651-685.         [ Links ]

[15] D. Sandri, Sur l'approximation numérique des écoulements quasi-newtoniens dont la viscosité suit la loi puissance ou la loi de Carreau, M2AN 27 (1993), 131-155.         [ Links ]

[16] R. Temam, Navier-Stokes equations, North-Holland, Amsterdam, (1977).         [ Links ]

[17] E. Zeidler, Nonlinear functional analysis and its applications II, Springer-Verlag, (1985).         [ Links ]