Applied Numerical Mathematics 181 (2022) 364–387

Contents lists available at ScienceDirect

Applied Numerical Mathematics

www.elsevier.com/locate/apnum

An adaptive iterative linearised finite element method for

implicitly constituted incompressible fluid flow problems and
its application to Bingham fluids
Pascal Heid ∗,1 , Endre Süli
Mathematical Institute, University of Oxford, Woodstock Road, Oxford OX2 6GG, UK

a r t i c l e i n f o a b s t r a c t

Article history: In this work, we further extend the theory in [C. Kreuzer, E. Süli, Adaptive finite element
Received 15 October 2021 approximation of steady flows of incompressible fluids with implicit power-law-like rheology,
Received in revised form 25 May 2022 ESAIM: Math. Model. Numer. Anal. 50 (5) (2016) 1333–1369] on the adaptive finite element
Accepted 17 June 2022
analysis of implicitly constituted incompressible fluid flow problems by taking into account
Available online 22 June 2022
the approximation of the nonlinear finite element solutions by an iterative solver. Thereby
Keywords: we show that the computable sequence generated by the modified algorithm still has a
Implicitly constituted incompressible fluid weak limit point, which is a solution of the given problem. Our abstract theory can be
flow problems applied to Bingham fluids, both with and without the convective term. The performance
Bingham fluids of the adaptive iterative linearised finite element algorithm in this context is illustrated by
Finite element methods two numerical experiments.
Kačanov scheme © 2022 The Author(s). Published by Elsevier B.V. on behalf of IMACS. This is an open
Zarantonello iteration
access article under the CC BY-NC-ND license
Adaptive algorithm
(http://creativecommons.org/licenses/by-nc-nd/4.0/).

1. Introduction

In this work, we consider implicitly constituted incompressible fluid flow problems, as introduced by Rajagopal in [37,
38]: Instead of demanding, as in the classical theory of continuum mechanics, that the Cauchy stress is an explicit function
of the symmetric part of the gradient of the velocity vector, one may allow an implicit relation of those two quantities. For
a rigorous mathematical analysis, including the proof of the existence of a weak solution for models of implicitly constituted
fluids we refer the reader to [9] and [10] for steady and unsteady flows, respectively. Regarding the numerical analysis of
such fluid flow models, very few results have been published so far. In [13] the authors proved the weak convergence of a
sequence of finite element approximations to a weak solution in the steady case. An a posteriori analysis of finite element
approximations of the model was carried out in [29]; in addition, those authors proved the weak convergence of an adaptive
finite element method. We further refer to [27] for a numerical investigation of several finite element discretisation methods
and linearisation schemes for Bingham and stress-power-law fluids. Just recently, in [19], the author proposed a semismooth
Newton method for the numerical approximation of a steady solution to implicit constitutive fluid flow problems. For the
numerical analysis of the unsteady case (not considered in this work) we refer to [16,39].

* Corresponding author.
E-mail addresses: pascal.heid@maths.ox.ac.uk (P. Heid), endre.suli@maths.ox.ac.uk (E. Süli).
1
PH acknowledges the financial support of the Swiss National Science Foundation (SNF), Project No. P2BEP2_191760.

https://doi.org/10.1016/j.apnum.2022.06.011
0168-9274/© 2022 The Author(s). Published by Elsevier B.V. on behalf of IMACS. This is an open access article under the CC BY-NC-ND license
(http://creativecommons.org/licenses/by-nc-nd/4.0/).
P. Heid and E. Süli Applied Numerical Mathematics 181 (2022) 364–387

The aim of this work is to extend the theory of [29] to adaptive iterative linearised finite element methods for implicitly
constituted fluid flow problems; i.e., in contrast to [29], we will take into account the approximation of a finite element
solution by a (linear) iteration scheme. Subsequently, it shall be shown that this abstract analysis can be applied in the
context of Bingham fluids.
The organisation of the paper is as follows: In Section 2, we present the preliminaries, including the formulation of
the problem, its finite element approximation, as well as an a posteriori error analysis. In Section 3, we will state the
adaptive iterative linearised finite element method (AILFEM), and prove the convergence of this algorithm. Since this result
can be shown by some minor modifications of the analysis in [29], we will omit the details of the proof, and rather
give a rough sketch in Appendix A and highlight where the main differences occur. Next, in Section 4, we will verify
the assumptions required for the convergence of AILFEM for Bingham flows. We will perform corresponding numerical
experiments in Section 5, and conclude our work with some closing remarks.

2. Preliminaries

In this section, we will introduce the model of a steady flow of an incompressible fluid in a bounded open Lipschitz do-
main  ⊂ Rd , d ∈ {2, 3}, with polyhedral boundary ∂ , which satisfies an implicit constitutive relation given by a maximal
monotone r-graph. Beforehand, let us introduce some basic notions concerning Lebesgue and Sobolev spaces.

2.1. Basic notations

For any measurable open set ω ⊆ , we denote by D (ω) := C 0∞ (ω; R) the space of smooth functions with compact
support in ω . Its dual space, i.e., the space of all bounded linear functionals on D (ω), is denoted by D  (ω). Similarly, we
denote by D (ω)d := C 0∞ (ω; Rd ) the space of vector-valued smooth functions with compact support in ω , and by D  (ω)d its
dual space.
Next, for any s ∈ [1, ∞) we denote by L s (ω) := L s (ω; R) the Lebesgue space of s-integrable functions with corresponding
 1/s
norm  f s,ω := ω | f (x)|s dx . Moreover, L ∞ (ω) = L ∞ (ω; R) denotes the  Lebesgue space of essentially bounded functions
with the norm  f ∞,ω := ess supx∈ω | f (x)|, and L 0s (ω) := { f ∈ L s (ω) : ω f dx = 0} denotes the set of functions (in the
 
corresponding Lebesgue space) with zero mean value. We note that, for s ∈ (1, ∞), L s (ω) and L 0s (ω) are the dual spaces of
L s (ω) and L 0s (ω), respectively, where s ∈ (1, ∞) is the Hölder conjugate of s, i.e., 1/s + 1/s = 1.
Likewise, for s ∈ [1, ∞], we denote by W 1,s (ω)d := W 1,s (ω; Rd ) the space of vector-valued Sobolev functions. Moreover,
1, s
for s ∈ [1, ∞], the space of vector-valued Sobolev functions with zero trace along the boundary is denoted by W 0 (ω)d
1, s
and is equipped with the norm f1,s,ω := ∇ fs,ω . Equivalently, for s ∈ [1, ∞), W 0 (ω)d is the closure in W 1,s (ω)d of
D(ω)d := C 0∞ (ω)d and its dual space, for any s ∈ (1, ∞), is denoted by W −1,s (ω)d . The norm in the dual space W −1,s (ω)d


is as usual given by

ϕ −1,s ,ω := sup ϕ, v ,
1,s
v∈ W
0
(ω)d
v1,s,ω =1

where ·, · signifies the duality pairing. In the sequel, for ω = , we omit the domain in the subscript of the norms; e.g.,
we write ·s := ·s, .
×d we denote by δ : κ the Frobenius inner-product, and by |κ | the Frobenius norm.
Finally, for δ, κ ∈ Rdsym

2.2. Problem formulation

As before, let  ⊂ Rd , d ∈ {2, 3}, be a bounded open Lipschitz domain with polyhedral boundary ∂ . For r ∈ (1, ∞), we
set

dr
2(d−r )
if r ≤ d3d
+2 ,
r̃ :=
r otherwise,

where r  ∈ (1, ∞) denotes the Hölder conjugate of r. Then, the implicitly constituted incompressible fluid flow problem
 
under consideration reads as follows: for f ∈ L r ()d find (u, p , S) ∈ W 0 ()d × L r̃0 () × L r ()d×d such that
1,r

div(u ⊗ u + pI − S) = f in D ()d , (1a)


div u = 0 in D (), (1b)
(Du(x), S(x)) ∈ A(x) for a.e. x ∈ ; (1c)

here, Du := 1 ×d := { ∈ Rd×d :
(∇ u + (∇ u)T ) ∈ Rdsym = T } signifies the symmetric velocity gradient. Moreover,
κ κ κ A:→
2
×d
Rdsym d×d
× Rsym is a maximal monotone r-graph, i.e., for almost every x ∈  the following properties are satisfied:

365
P. Heid and E. Süli Applied Numerical Mathematics 181 (2022) 364–387

(A1) (0, 0) ∈ A(x);

(A2) For all (κ 1 , δ 1 ), (κ 2 , δ 2 ) ∈ A(x),

(δ 1 − δ 2 ) : (κ 1 − κ 2 ) ≥ 0;
×d × Rd×d and
(A3) If (κ , δ) ∈ Rdsym sym

(δ − δ) : (κ − κ ) ≥ 0 for all (κ , δ) ∈ A(x),

then (κ , δ) ∈ A(x);
(A4) There exist a non-negative function m ∈ L 1 () and a constant c > 0, such that for all (κ , δ) ∈ A(x) we have that

δ : κ ≥ −m(x) + c (|κ |r + |δ|r );
(A5) The set-valued mapping A :  → Rdsym ×d × Rd×d is measurable, i.e., for any closed sets C , C ⊂ Rd×d we have that
sym 1 2 sym
({x ∈  : A(x) ∩ (C1 × C2 ) = ∅}) is a Lebesgue measurable subset of .

The existence of a (not necessarily unique) solution to (1) for r > d2d
+2 was first established in [9]. The proof is based on
×d → Rd×d of the graph A; i.e., for all κ ∈ Rd×d we have that (κ , S (x, κ )) ∈
a (monotone) measurable selection S :  × Rdsym sym sym
×d → Rd×d was approximated by a sequence
A(x) for almost every x ∈ . In turn, the measurable selection S :  × Rdsym sym
d×d d×d
of strictly monotone mappings S :  × Rsym → Rsym , n = 0, 1, 2, . . . , obtained by a mollification of S . In this work,
n

following [29], we allow for more general graph approximations.

×d → Rd×d such that

Assumption 2.1. For any n ∈ N there exists a mapping Sn :  × Rdsym sym

×d is measurable for all κ ∈ Rd×d ;

• Sn (·, κ ) :  → Rdsym sym
d×d ×d is continuous for almost every x ∈ ;
• S (x, ·) : Rsym → Rdsym
n

×d , we have that
• Sn is strictly monotone in the sense that, for all κ 1 = κ 2 ∈ Rdsym

(Sn (x, κ 1 ) − Sn (x, κ 2 )) : (κ 1 − κ 2 ) > 0 for almost every x ∈ ;


• There exist constants c̃ 1 , c̃ 2 > 0 and non-negative functions m̃ ∈ L 1 () and k̃ ∈ L r () such that, uniformly in n ∈ N ,

|Sn (x, κ )| ≤ c̃ 1 |κ |r −1 + k̃(x), (2)

n r
S (x, κ ) : κ ≥ c̃ 2 |κ | − m̃(x) (3)

for all κ∈ ×d
Rdsym and almost every x ∈ .

Of course, we also need that Sn approximates, in a certain sense, the measurable selection S of A; this will be made

precise in Section 3. Then, the regularised counterpart of problem (1) is given as follows: for f ∈ L r ()d find (un , pn ) ∈
1,r
W 0 ()d × L r̃0 () such that

div(un ⊗ un + pn I − Sn (un )) = f in D ()d ,

div un = 0 in D ().
We note that in the problem formulation above and in the following, the explicit dependence of Sn on x ∈  will be
suppressed.

2.3. Finite element spaces

In this work, we consider a sequence {T N } N of shape-regular conforming triangulations of , such that T N +1 is obtained
by a refinement of T N . For any m ∈ N let us denote by Pm the space of polynomials of degree at most m. Then, the
corresponding conforming finite element spaces are given by

V (T N ) := {V ∈ W 01,2 ()d : V| K ∈ PV ( K ) for all K ∈ T N },

Q(T N ) := { Q ∈ L ∞ () : Q | K ∈ PQ ( K ) for all K ∈ T N },
where PV and PQ are spaces of polynomials such that P1d ⊆ PV ⊆ Pid and P0 ⊆ PQ ⊆ P j  for some i  ≥ j  ≥ 0; here,

Pmd = Pm × . . . × Pm , m ∈ N,
 
d times

366