Finite element approximation of steady flows of generalized Newtonian fluids with concentration-dependent power-law index
Seungchan Ko Endre Suli









Finite element approximation of steady flows of generalized
              Newtonian fluids with concentration-dependent power-law index
    Seungchan KoMathematical Institute, University of Oxford, Andrew Wiles Building, Woodstock Road, Oxford OX2 6GG, UK. Email: seungchan.ko@maths.ox.ac.uk
 and  Endre SüliMathematical Institute, University of Oxford, Andrew Wiles Building, Woodstock Road, Oxford OX2 6GG, UK. Email: endre.suli@maths.ox.ac.uk
     
=====================================================================================================================================================================================================================================================================================================================



 



We consider a system of nonlinear partial differential equations describing the motion of an incompressible chemically reacting
generalized Newtonian fluid in three space dimensions. The governing system consists of a steady convection-diffusion equation
for the concentration and a  generalized steady power-law-type fluid flow model for the velocity and the pressure, where the
viscosity depends on both the shear-rate and the concentration through a concentration-dependent power-law index.
The aim of the paper is to perform a mathematical analysis of a finite element approximation of this model. We formulate a
regularization of the model by introducing an additional term in the conservation-of-momentum equation and construct a
finite element approximation of the regularized system. We show the convergence of the finite element method to a weak solution
of the regularized model and prove that weak solutions of the regularized problem converge to a weak solution of the
original problem.


Keywords: Non-Newtonian fluid, variable exponent, synovial fluid, finite element method



AMS Classification: 65N30, 74S05, 76A05


Introduction
We are interested in developing a convergence theory for finite element approximations of a system of nonlinear partial differential
equations (PDEs) modelling the rheological response of the synovial fluid. The synovial fluid is a biological fluid found in the
cavities of movable joints and is composed of ultrafiltrated blood, called hyaluronan. Laboratory experiments have shown that
the viscosity of the fluid depends on the concentration of hyaluronan, as well as on the shear-rate. In particular, it was observed
in steady shear experiments that the concentration of the hyaluronan is not just a scaling factor of the viscosity (understood as
ν(c,|Du|)=f(c) ν̃(|Du|)) but it has an influence on the degree of shear-thinning. Therefore, a new mathematical model
of the rheological response of the synovial fluid was proposed in <cit.>. There, the authors considered a power-law-type model
for the velocity and the pressure, where the power-law index depends on the concentration, corresponding to the fact that the
concentration affects the level of shear-thinning. To close the system, a generalized convection-diffusion equation was assumed to
be satisfied by the concentration. For a detailed rheological background we refer to <cit.>.

Based on the description above, we consider the following system of PDEs:


div    =0          in  Ω,

div (⊗)-div (c,)   =-∇p+f          in  Ω,

div (c)-div _c(c,∇c, )   =0          in  Ω,


where Ω⊂^d is a bounded open Lipschitz domain. In the above system of PDEs,
:Ω→^d, p:Ω→, c:Ω→_≥0 denote the velocity,
pressure and concentration fields, respectively, f:Ω→^d is a given external force, and  denotes
the symmetric velocity gradient, i.e. =1/2(∇+(∇)^T). To complete the model, we impose the following
Dirichlet boundary conditions:


    =0,     c=c_d    on  ∂Ω,


where c_d∈ W^1,s(Ω) for some s>d. By Sobolev embedding, c_d is continuous up to the boundary, and we can therefore define


    c^-min_x∈Ωc_d     and     c^+max_x∈Ωc_d.



We further assume that the extra stress tensor :_≥ 0×^d×
d_ sym→^d× d_ sym is a continuous mapping
satisfying the following growth, strict monotonicity and coercivity
conditions, respectively: there exist positive constants C_1, C_2 and C_3 such that


    |(ξ,B)|≤ C_1(|B|^r(ξ)-1+1),



    ((ξ,B_1)-(ξ,B_2))·(B_1-B_2)>0   for  B_1≠B_2,



    (ξ,B)·B≥ C_2(|B|^r(ξ)+||^r'(ξ))-C_3,


where r:_≥0→_≥0 is a Hölder-continuous function satisfying 1<r^-≤ r(ξ)≤ r^+<∞ and r'(ξ)
is defined as its Hölder conjugate, r(ξ)/r(ξ)-1. We further assume that the concentration flux vector
_c(ξ,g,B):_≥ 0×^d×^d× d_sym→^d is a
continuous mapping, which is linear with respect to g, and it additionally satisfies the following growth and coercivity
conditions: there exist positive constants C_4 and C_5 such that


    |_c(ξ,g,B)|   ≤ C_4|g|,
    _c(ξ,g,B)·g   ≥ C_5|g|^2.



As we have discussed above, the prototypical examples we have in mind are the following:


    (c,)=ν(c,||),    _c(c,∇ c,)=K(c,||)∇ c,


where the viscosity ν(c,||) is of the following form:

    ν(c,||)∼ν_0(κ_1+κ_2||^2)^r(c)-2/2

and ν_0,κ_1,κ_2 are positive constants.

The rigorous mathematical analysis of the existence of global weak solutions to a PDE system, consisting of the generalized
Navier–Stokes equations, with a concentration-dependent viscosity coefficient, coupled to a convection-diffusion equation,
was initiated in <cit.>. There, however, the power-law index was fixed and the concentration was only a scaling factor of the
viscosity; the authors considered the evolutionary model and established the long-time existence of large-data global weak solutions.
Concerning the model (<ref>)–(<ref>) where the power-law index is concentration-dependent, the mathematical analysis was
initiated in <cit.>. The authors established there the existence of weak solutions, provided that r^->3d/d+2,
by using generalized monotone operator theory. In <cit.>, with the help of a Lipschitz-truncation technique, the existence
of weak solutions with r^->d/2 was proved and the Hölder continuity of the concentration was shown by using De Giorgi's method.

In <cit.>, the convergence of a finite element approximation to the system (<ref>)–(<ref>) was shown, using a discrete
De Giorgi regularity result. Because of the absence of a discrete De Giorgi regularity result in three space dimensions, the analysis in
<cit.> was restricted to the case of two space dimensions. In this paper, we extend the analysis developed in <cit.> to
three space dimensions, and we formulate an analogous convergence result for a finite element method in a three-dimensional domain.
The main idea here is to use a different numerical approximation scheme from the one in <cit.>, resulting in a different
limiting process. To this end, we consider different meshes for the conservation of linear momentum equation and the concentration
equation. The resulting numerical method can be viewed as a two-level Galerkin approximation. This enables us to
separate the passages to the limits with respect to the discretization parameters in the two equations, thus avoiding the need for
a discrete De Giorgi regularity result in three space dimensions.

As a first step, in Section 2 we introduce the necessary notational conventions and auxiliary results, which will be used throughout the paper.
In Section 3, we define a regularized problem, which enables us to enlarge the range of the power-law index so as to be able to cover the
practically relevant range of values of this index. In Sections 4 and 5, we construct a two-level Galerkin finite element approximation to
the regularized problem and perform a convergence analysis of the numerical method. Finally, in Section 6, we shall prove that weak
solutions of the regularized problem converge to a weak solution of the original problem when we pass to the limit with the regularization
parameter.




Notation and auxiliary results
In this section, we shall introduce certain function spaces and auxiliary results that will be used throughout the paper. Let 𝒫 be
the set of all measurable functions r:Ω→[1,∞]; we shall call the function r∈𝒫(Ω) a variable
exponent. We define r^- ess inf_x∈Ωr(x), r^+ ess sup_x∈Ωr(x) and we only consider the case


    1<r^-≤ r^+<∞.



Since we are considering a power-law index depending on the concentration, we need to work with Lebesgue and Sobolev spaces with
variable exponents. To be specific, we introduce the following variable-exponent Lebesgue spaces, equipped with the corresponding
Luxembourg norms:


    L^r(·)(Ω)   {u∈ L^1_loc(Ω):∫_Ω|u(x)|^r(x)<∞},
    u_L^r(·)(Ω)=u_r(·)   inf{λ>0:∫_Ω|u(x)/λ|^r(x)≤1}.


Similarly, we introduce the following generalized Sobolev spaces:


    W^1,r(·)(Ω)   {u∈ W^1,1(Ω)∩ L^r(·)(Ω):|∇ u|∈ L^r(·)},
    u_W^1,r(·)(Ω)=u_1,r(·)   inf{λ>0:∫_Ω[|u(x)/λ|^r(x)+|∇ u(x)/λ|^r(x)]≤1}.


It is easy to show that all of the above spaces are Banach spaces, and because of (<ref>), they are all separable and
reflexive; see <cit.>.

Furthermore, we introduce certain function spaces that are frequently used in PDE models of incompressible fluids. Henceforth,
X(Ω)^d will denote the space of d-component vector-valued functions with components from X(Ω). We also define
the space of tensor-valued functions X(Ω)^d× d. Finally, we define the following spaces:


    W^1,r(·)_0(Ω)^d   {∈ W^1,r(·)(Ω)^d:=0  on  ∂Ω},
    
    W^1,r(·)_0,div(Ω)^d   {∈ W^1,r(·)_0(Ω)^d:div =0   in Ω},
    
    L^r(·)_0(Ω)   {f∈ L^r(·)(Ω):∫_Ωf(x)=0}.



Throughout the paper, we shall denote the duality pairing between f∈ X and g∈ X^* by  ⟨ g,f⟩, and for two
vectors a and b, a·b denotes their scalar product; similarly, for
two tensors 𝐀 and 𝐁, 𝐀·𝐁 signifies their scalar product. Also, for any Lebesgue-measurable set Q⊂^d, |Q| denotes the standard Lebesgue measure of the set Q.

Next we introduce the necessary technical tools. First we define the subset 𝒫^log(Ω)⊂𝒫(Ω):
it will denote the set of all log-Hölder-continuous functions
defined on Ω, that is the set of all functions r ∈𝒫(Ω) satisfying


    |r(x)-r(y)|≤C_log(r)/-log|x-y|    ∀  x,y∈Ω:0<|x-y|≤1/2.


It is obvious that classical Hölder-continuous functions on Ω automatically belong to this class.

Next we state the following lemma, which summarizes some inequalities involving variable-exponent norms.
For proofs, see <cit.>, which is an extensive source of information concerning variable-exponent spaces.



Let Ω⊂^d be a bounded open Lipschitz domain and let r∈𝒫^log(Ω) satisfy (<ref>).
Then, the following inequalities hold:




  * Hölder's inequality, i.e.,


    fg_s(·)≤2f_r(·)g_q(·),   with   r,q,s∈𝒫(Ω),
       1/s(x)=1/r(x)+1/q(x),    x ∈Ω.




  * Poincaré's inequality, i.e.,


    u_r(·)≤ C(d,C_log(r)) diam(Ω)∇ u_r(·)    ∀  u∈ W^1,r(·)_0(Ω).




  * Korn's inequality, i.e.,


    ∇_r(·)≤ C(Ω, C_log(r))_r(·)    ∀ ∈ W^1,r(·)_0(Ω)^d,


where C_log(r) is the constant appearing in the definition of the class of log-Hölder-continuous functions.





Another important auxiliary result is the existence of the Bogovskiĭ operator in the variable-exponent setting.

Let Ω⊂^d be a bounded open Lipschitz domain and suppose that r∈𝒫^log(Ω) with 1<r^-≤ r^+<∞.
Then, there exists a bounded linear operator ℬ:L^r(·)_0(Ω)→ W^1,r(·)_0(Ω)^d such that for
all f∈ L^r(·)_0(Ω) we have


    div (ℬf)   =f,
    ℬf_1,r(·)   ≤ Cf_r(·),


where C depends on Ω, r^-, r^+, and C_log(r).



Let us now state the inf-sup condition, which has a crucial role in the mathematical analysis of incompressible fluid flow problems.


For any s, s'∈(1,∞), with 1/s+1/s'=1, there exists a positive constant α_s>0 such that

    α_sq_s'≤sup_0≠∈ W^1,s_0(Ω)^d⟨div ,q⟩/_1,s    ∀ 
    q∈ L^s'_0(Ω).


This is a direct consequence of the existence of the Bogovskiĭ operator in spaces with fixed exponent, which is a special
case of Theorem <ref>;  see <cit.> for additional details.

Furthermore, we can prove the following inf-sup condition in spaces with variable-exponent norms, which will play an important role
in the subsequent analysis.


Let Ω⊂^d be a bounded open Lipschitz domain and let r∈𝒫^log(Ω) with 1<r^-≤ r^+<∞.
Then, there exists a constant α_r>0 such that

    α_rq_r'(·)≤sup_0≠∈ W^1,r(·)_0(Ω)^d⟨div ,q⟩/_1,r(·)    ∀  q∈ L^r'(·)_0(Ω).



Proposition <ref> is a direct consequence of Theorem <ref> and the norm-conjugate formula stated in the following lemma.


Let r∈𝒫^log(Ω) be a variable exponent with 1<r^-≤ r^+<∞; then we have

    1/2f_r(·)≤sup_g∈ L^r'(·)(Ω),  g_r'(·)≤1∫_Ω|f||g|,

for all measurable functions f ∈  L^r(·)(Ω).


Finally, we recall the following well-known result due to De Giorgi and Nash <cit.>; see also <cit.> for its
application to the system of partial differential equations considered in the present paper.


Let Ω⊂^d be a Lipschitz domain and let s>d be fixed. Suppose that K∈ L^∞(Ω)^d× d
is uniformly elliptic with ellipticity constant λ>0. Then, there exists an α∈(0,1) such that, for any
f∈ L^s(Ω)^d, g∈ L^ds/d+s(Ω) and any c_d∈ W^1,s(Ω), there exists a unique
c∈ W^1,2(Ω) such that c-c_d∈ W^1,2_0(Ω)∩ C^0,α(Ω) and


    ∫_ΩK∇ c·∇φ=∫_Ωf·∇φ+∫_Ωgφ    ∀ φ∈ W^1,2_0(Ω);


furthermore, the following uniform bound holds:

    c_W^1,2∩ C^0,α≤ C(Ω,λ,s,K_∞,f_s,g_ds/d+s,
    c_d_1,s).





Using these notations, the weak formulation of the problem (<ref>)–(<ref>) is as follows.


Problem (Q). For f∈ (W^1,r^-_0(Ω)^d)^*, c_d∈ W^1,s(Ω), s>d, and a Hölder-continuous
function r, with 1<r^-≤ r(c)≤ r^+<∞ for all c∈[c^-,c^+], find (c-c_d)∈ W^1,2_0(Ω)∩
C^0,α(Ω), for some α∈ (0,1), ∈ W^1,r(c)_0(Ω)^d, p∈ L^r'(c)_0(Ω) such that


∫_Ω(c,)·∇ψ-(⊗)·∇ψ-⟨div 
ψ,p⟩   =⟨f,ψ⟩          ∀ ψ∈W^1,∞_0(Ω)^d,

∫_Ωq div    =0          ∀  q∈L^r'(c)_0(Ω),

∫_Ω_c(c,∇c,)·∇φ-c·∇φ   =0          ∀ φ∈W^1,2_0(Ω).



Thanks to Proposition <ref>, we can restate Problem (Q) in the following (equivalent) divergence-free setting.


Problem (P). For f∈ (W^1,r^-_0(Ω)^d)^*, c_d∈ W^1,s(Ω), s>d, and a Hölder-continuous
function r, with 1<r^-≤ r(c)≤ r^+<∞ for all c∈[c^-,c^+], find (c-c_d)∈ C^0,α(Ω)
∩ W^1,2_0(Ω), ∈ W^1,r(c)_0,div(Ω)^d, such that


    ∫_Ω(c,)·∇ψ-(⊗)·∇ψ   =⟨f,ψ⟩      ∀ ψ∈ W^1,∞_0,div(Ω)^d,
    ∫_Ω_c(c,∇ c,)·∇φ-c·∇φ   =0      ∀ φ∈ W^1,2_0(Ω).



From now on, for simplicity, we shall restrict ourselves to the case of d=3. Our results can be however easily extended to
the case of any d ≥ 2. We note in passing that since no uniqueness result is currently known for weak solutions of the
problem under consideration, we can only prove that a subsequence of the sequence of discrete solutions converges to a weak
solution of the problem.




Regularization of the problem
Before constructing the approximation of problem (Q) we shall formulate a regularized problem; it will then be the
regularized problem that will be approximated by a finite element method. We shall show that the sequence of finite element
approximations converges to a weak solution of the regularized problem, and that solutions of the regularized problem, in turn,
converge to a weak solution of problem (Q). The reason for proceeding in this way is that direct approximation of problem (Q),
which bypasses the use of the regularized problem, necessitates the imposition of an unnaturally strong condition on the variable exponent
r in the convergence analysis of the finite element method; the procedure that we describe below does not suffer from this
shortcoming.

Motivated by <cit.>,
we shall utilize the following regularized problem, involving the regularization parameter k∈ℕ.
We choose a sufficiently large t>0, such that r^->3/2>t/t-2. Then we seek a weak solution
(,p,c)(^k,p^k,c^k) to


    div    =0          in  Ω,
    div (⊗)-div (c,)+1/k||^t-2   =
    -∇ p+f          in  Ω,
    div (c)-div _c(c,∇ c, )   =0          in  Ω,


Therefore, we consider the following regularized weak formulation.

Problem (Q*). For f∈ (W^1,r^-_0(Ω)^3)^*, c_d∈ W^1,s(Ω), s>3, and a Hölder-continuous
function r, with 1<r^-≤ r(c)≤ r^+<∞ for all c∈[c^-,c^+], and r^->3/2>t/t-2, t>2, find
(c-c_d)(c^k-c_d)∈ W^1,2_0(Ω)
∩ C^0,α(Ω), for some α∈ (0,1), ^k∈ W^1,r(c)_0(Ω)^3, p
p^k∈ L^r'(c)_0(Ω) such that


    ∫_Ω(c,)·∇ψ-(⊗)·∇ψ+
    1/k||^t-2·ψ-⟨div ψ,p⟩   =
    ⟨f,ψ⟩        ∀ ψ∈ W^1,∞_0(Ω)^3,
    ∫_Ωq div    =0        ∀  q∈ L^r'(c)_0(Ω),
    ∫_Ω_c(c,∇ c,)·∇φ-c·∇φ   =0        ∀ φ∈ W^1,2_0(Ω).



Again, by using Proposition <ref>, we can restate Problem (Q*) in the following (equivalent) divergence-free setting:

Problem (P*). For f∈ (W^1,r^-_0(Ω)^3)^*, c_d∈ W^1,s(Ω), s>3, and Hölder-continuous function
r, with 1<r^-≤ r(c)≤ r^+<∞ for all c∈[c^-,c^+],  and r^->3/2>t/t-2, t>2, find (c-c_d)(c^k-c_d)∈ C^0,α(Ω)∩
W^1,2_0(Ω), ^k∈ W^1,r(c)_0,div(Ω)^3, such that


    ∫_Ω(c,)·∇ψ-(⊗)·∇ψ+1/k||^t-2·ψ   =⟨f,ψ⟩      ∀ ψ∈ W^1,∞_0,div(Ω)^3,
    ∫_Ω_c(c,∇ c,)·∇φ-c·∇φ   =0      ∀ φ∈ W^1,2_0(Ω).



We shall formulate the finite element approximation of the regularized problem Problem (Q*) in a three-dimensional domain;
the convergence analysis of the method
is presented in Section 4 and Section 5. In Section 6, we will prove that a sequence of weak solution triples
{(^k,p^k,c^k)}_k ≥ 1
of the regularized problem converges to a weak solution triple (,p,c) of Problem (Q).
The latter result is recorded in our next theorem.



Suppose that Ω⊂^3 is a convex polyhedral domain and c_d∈ W^1,s(Ω) for some s>3.
Let us further assume that r:_≥0→_≥0 is a Hölder-continuous function with r^->3/2>t/t-2, t>2,
and suppose that f∈ (W^1,r^-_0(Ω)^3)^*.
Let (^k,p^k,c^k) be a weak solution of the regularized problem (<ref>)–(<ref>).
Then, as k→∞, (a subsequence, not indicated, of) the sequence {(^k,p^k,c^k)}_k ≥ 1
converges to (,p,c) in the following sense:


    ^k    ⇀      weakly  in  W^1,r^-_0,div(Ω)^3, 
    
    c^k    ⇀ c       weakly  in  W^1,2(Ω),
    
    c^k    → c       strongly  in  C^0,α(Ω)    for   some   α∈(0,1),
    
    p^k    ⇀ p       weakly  in  L^j'(Ω)    ∀  j>max{r^+,2}.

Furthermore, (,p,c) is a weak solution of the problem Problem (Q*) stated in (<ref>)–(<ref>).








Finite element approximation
Finite element spaces
Let {_n}, {_m}  be families of shape-regular partitions of Ω such that the following properties hold:


  * Affine equivalence: For each element E∈_n (or E∈_m) , there exists an invertible affine mapping

    F_E:E→Ê,

where Ê is the standard reference 3-simplex in ^3.

  * Shape-regularity: For any element E∈_n (or E∈_m), the ratio of diam E to
the radius of the inscribed ball is bounded below uniformly by a positive constant, with respect to all _n (or _m) and
n∈ℕ (or m ∈ℕ).




For given partitions _n and _m, the finite element spaces are defined by


    ^n   =(_n){∈ C(Ω)^3:_|E∘F^-1_E∈ℙ̂_,E∈_n  and  _|∂Ω=0},
    ^n   =(_n){Q∈ L^∞(Ω):Q_|E∘F^-1_E∈ℙ̂_,E∈_n},
    ^m   =(_m){Z∈ C(Ω):Z_|E∘F^-1_E∈ℙ̂_,E∈_m  and  Z_|∂Ω=0},


where ℙ̂_⊂ W^1,∞(Ê)^3, ℙ̂_⊂ L^∞(Ê) and ℙ̂_⊂ W^1,∞(Ê) are finite-dimensional linear subspaces.

We assume that ^n and ^m have finite and locally supported bases; for example,
for each n∈ℕ and m∈ℕ, there exists an N_n∈ℕ and an N_m∈ℕ such that


    ^n=span{^n_1,…,^n_N_n},


    ^m=span{Z^m_1,…,Z^m_N_m},


and for each basis function ^n_i, Z^m_j, we have that if there exists an E∈_n (respectively, _m), with ^n_i≠0 (respectively, Z^m_j≠0) on E, then


    supp ^n_i⊂⋃{E'∈_n:E'∩ E≠∅} S_E.


    supp Z^m_j⊂⋃{E'∈_m:E'∩ E≠∅} T_E.


For the pressure space ^n, we assume that ^n has a basis consisting of discontinuous piecewise polynomials; i.e., for each n∈ℕ, there exists an Ñ_n∈ℕ such that


    ^n=span{Q^n_1,…,Q^n_Ñ_n}


and for each basis function Q^n_i, we have that


    supp Q^n_i=E    for   some   E∈_n.


We assume further that ^n contains continuous piecewise linear functions and ^n contains piecewise constant functions.

Using the assumed shape-regularity we can easily verify that


    ∃ X∈ℕ:|S_E|≤ X|E|   for   all   E∈_n,


    ∃ Y∈ℕ:|T_E|≤ Y|E|   for   all   E∈_m,


where X is independent of n and Y is independent of m. We denote by g_E the diameter of E∈_n and
by h_E the diameter of E∈_m.

We also introduce the subspace ^n_div of discretely divergence-free functions. More precisely, we define

    ^n_div{∈^n:⟨div ,Q⟩=0   ∀  Q∈^n},

and the subspace of ^n consisting of vanishing integral mean-value approximations:

    ^n_0{Q∈^n:∫_ΩQ=0}.


Throughout this paper, we assume that the finite element spaces introduced above have the following minimal approximation properties.



Assumption 1 (Approximability) For all s∈[1,∞),

    inf_∈^n-_1,s   →0        ∀ ∈ W^1,s_0(Ω)^3  as  n→∞, 
    inf_Q∈^nq-Q_s   →0         ∀  q∈ L^s(Ω)  as  n→∞,
    inf_Z∈^mz-Z_1,s   →0        ∀  z∈ W^1,s_0(Ω)  as  m→∞.

For this, a necessary condition is that the maximal mesh size vanishes, i.e., that max_E∈_ng_E→0 as n→∞ and max_E∈_mh_E→0 as m→∞.



Assumption 2 (Existence of a projection operator Π^n_div) For each n∈ℕ, there exists a linear projection operator Π^n_div:W^1,1_0(Ω)^3→^n such that:


  * Π^n_div preserves the divergence structure in the dual of the discrete pressure space; in other words, for any ∈ W^1,1_0(Ω)^3, we have

    ⟨div ,Q⟩=⟨div Π^n_div,Q⟩    ∀  Q∈^n.


  * Π^n_div is locally W^1,1-stable, i.e., there exists a constant c_1>0, independent of n, such that

    _E|Π^n_div|+g_E|∇Π^n_div|≤ c_1_S_E||+g_E|∇|    ∀ ∈ W^1,1_0(Ω)^3   and   ∀  E∈_n.






Note that the local W^1,1(Ω)^3-stability of Π^n_div implies its local and global W^1,s(Ω)^3-stability for s∈[1,∞]. In other words, for any s∈[1,∞] we have

    Π^n_div_1,s≤ c_s_1,s    ∀ ∈ W^1,s_0(Ω)^3,

with a constant c_s>0 independent of n>0.

 Note further that the approximability (Assumption 1) and inequality (<ref>) imply the convergence of Π^n_div to . In fact, 

    -Π^n_div_1,s→0    ∀ ∈ W^1,s_0(Ω)^3  as  n→∞,   ∀  s ∈ [1,∞).




Assumption 3 (Existence of a projection operator Π^n_) For each n∈ℕ,
there exists a linear projection operator Π^n_:L^1(Ω)→^n such that Π^n_ is locally L^1-stable; i.e.,
there exists a constant c_2>0, independent of n, such that

    _E|Π^n_q|≤ c_2_S_E|q|

for all q∈ L^1(Ω) and all E∈_n.



Again, we have the following global stability and convergence property:

    Π^n_q_s'≤ c_s'q_s'    ∀  q∈ L^s'(Ω),   ∀  s'∈(1,∞),

and

    q-Π^n_q_s'→0,    as   n→∞   for   all   q∈ L^s'(Ω)   and   s'∈(1,∞).



According to <cit.>, the following pairs of velocity-pressure finite element spaces satisfy Assumptions 1, 2 and 3, for example:


  * The conforming Crouzeix–Raviart Stokes element, i.e., continuous piecewise quadratic plus cubic bubble velocity and discontinuous piecewise linear pressure approximation (compare e.g. with <cit.>);

  * The space of continuous piecewise quadratic polynomials for the velocity and piecewise constant pressure approximation; see, <cit.>.




Our final assumption is the existence of a projection operator for the concentration space.



Assumption 4 (Existence of a projection operator Π^m_) For each m∈ℕ, there exists a linear projection operator Π^m_:W^1,1_0(Ω)→^m such that

    _E|Π^m_z|+h_E|∇Π^m_z|≤ c_3_T_E|z|+h_E|∇ z|    ∀  z∈ W^1,1_0(Ω)   and   ∀  E∈_m,

where c_3 does not depend on m.

Similarly as above, the projection operator Π^m_ is globally W^1,s-stable for s∈[1,∞], and thus, by approximability,

    Π^m_z-z_1,s→0     ∀  z∈ W^1,s_0(Ω),  ∀  s ∈ [1,∞).


Finally, we introduce a discrete inf-sup condition, which holds in our finite element setting.
It is a direct consequence of (<ref>) and the existence of Π^n_div; see
<cit.> for further details.


For s, s'∈(1,∞) satisfying 1/s+1/s'=1, there exists a positive constant β_r>0, which is independent of n, such that

    β_rQ_s'≤sup_0≠∈^n⟨div ,Q⟩/_1,s    ∀  Q∈^n_0   and   ∀ n∈ℕ.






The finite element approximation
In this section, we shall construct the finite element approximation of the problem (<ref>)–(<ref>).
An important property of the incompressible Navier–Stokes equations is that the convective term in the momentum equation
is skew-symmetric; this is a consequence of the velocity field  being divergence-free.
However, in the discretized problem, we might lose the skew-symmetry because we are considering only discretely divergence-free
finite element functions from the finite element space for the velocity.
Thus we need to modify the finite element approximation of the convective term in order to ensure that the skew-symmetry
is preserved under discretization. We therefore define the following modified convective terms:


    B_u[,,]   1/2∫_Ω((⊗)·∇-(⊗)·∇),
    
    B_c[b,,z]   1/2∫_Ω(z·∇ b-b·∇ z),


for all ,,∈ W^1,∞_0(Ω)^3, b,z∈ W^1,∞(Ω). These trilinear forms then
coincide with the corresponding trilinear forms appearing in the weak formulations of the momentum equation and the
concentration equation, provided that we are considering
pointwise divergence-free velocity fields. Furthermore, thanks to their skew symmetry, these two trilinear forms
now also vanish for discretely
divergence-free functions when = and b=z, respectively. Explicitly, we have


    B_u[,,]   =0     and     B_c[z,,z]=0      ∀ ∈ W^1,∞_0(Ω)^3,   z∈ W^1,∞(Ω),
    
    B_u[,,]   =-∫_Ω(⊗)·∇      ∀ ,,∈ W^1,∞_0,div(Ω)^3,
    
    B_c[b,,z]   =-∫_Ωb·∇ z      ∀ ∈ W^1,∞_0,div(Ω)^3,   b,z∈ W^1,∞(Ω).


Moreover, the trilinear form B_u[·,·,·] is bounded. Indeed, if ,,∈ W^1,∞_0(Ω)^3,
then, by Hölder's inequality,


    ∫_Ω(⊗)·∇≤_2(r^-)'_2(r^-)'_1,r^-,


and


    ∫_Ω(⊗)·∇≤_2(r^-)'_2(r^-)'_1,r^-.


Therefore, we obtain the bound


    |B_u[,,]|≤_2(r^-)'_2(r^-)'_1,r^-+_2(r^-)'_1,r^-_2(r^-)'.



Now, for each n,m∈ℕ, we call a triple (^n,m,P^n,m,C^n,m)∈^n×^n_0×(^m+c_d) a discrete
solution to the Galerkin approximation if it satisfies


∫_Ω(C^n,m,^n,m)·+1/k|^n,m|^t-2^n,m·   +B_u[^n,m,^n,m,]

-⟨div ,P^n,m⟩   =⟨f,⟩      ∀ ∈^n,

∫_ΩQ div ^n,m   =0      ∀  Q∈^n,

∫_Ω_c(C^n,m,∇C^n,m,^n,m)·∇Z+B_c[C^n,m,^n,m,Z]   =0      ∀  Z∈^m,


where c_d∈ W^1,s(Ω) with s>3 and f∈(W^1,r^-_0(Ω)^3)^*.

If we restrict the test functions  to ^n_div, then the above problem is transformed to the following:
find (^n,m,C^n,m)∈^n_div×(^m+c_d) satisfying


    ∫_Ω(C^n,m,^n,m)·+1/k|^n,m|^t-2^n,m·
    +B_u[^n,m,^n,m,]   =⟨f,⟩      ∀ ∈^n_div
    ∫_Ω_c(C^n,m,∇ C^n,m,^n,m)·∇ Z+B_c[C^n,m,^n,m,Z]   =0      ∀  Z∈^m.



If 3/2<r^-, the existence of the discrete solution pair (^n,m, C^n,m)∈^n_div×(^m+c_d)
follows from a fixed point argument combined with an iteration scheme. Let us briefly summarize the proof of the existence
of the pair (^n,m, C^n,m)∈^n_div×(^m+c_d). Let {_i}^N_n_i=1 be a basis of
^n_div⊂ W^1,∞_0(Ω)^3 such that ∫_Ω_i·_j=δ_ij and let
{z_j}^N_m_j=1 be a basis of ^m⊂ W^1,2_0(Ω) such that ∫_Ωz_iz_j=δ_ij.
Then, for fixed n,m∈ℕ, we define the Galerkin approximations.


    ^n,m∑^N_n_i=1α_i^n,m_i,     C^n,m∑^N_m_i=1β_i^n,mz_i+c_d,


which satisfy (<ref>)–(<ref>).

 First we define C^n,m_1 c_d∈^m+c_d. Then, for any ℓ∈ℕ, we define ^n,m_ℓ∈^n_div
 as a solution of the finite-dimensional problem
 

    ∫_Ω(C^n,m_ℓ,^n,m_ℓ)·+1/k|^n,m_ℓ|^t-2^n,m·+B_u[^n,m_ℓ,^n,m_ℓ,]=⟨f,⟩    ∀ ∈^n_div,


and C^n,m_ℓ∈^m+c_d as a solution of the finite-dimensional problem


    ∫_Ω_c(C^n,m_ℓ,∇ C^n,m_ℓ,^n,m_ℓ-1)·∇ Z+B_c[C^n,m_ℓ,^n,m_ℓ-1,Z]=0    ∀  Z∈^m.


The existence of the functions  ^n,m_ℓ∈^n_div and C^n,m_ℓ∈^m+c_d is easily shown by means of Brouwer's fixed point theorem. Furthermore, for each n,m∈ℕ, the sequences of functions {^n,m_ℓ}^∞_ℓ=1 and {C^n,m_ℓ}^∞_ℓ=1 satisfy the following uniform bounds:


    ^n,m_ℓ_1,r^-+^n,m_ℓ_t≤ C_1,    ∇ C^n,m_ℓ_2≤ C_2,


where C_1 and C_2 are positive constants, independent of ℓ.
Thus, by the Bolzano–Weierstrass theorem we deduce the existence of
limits ^n,m∈^n_div and C^n,m∈^m+c_d for ^n,m_ℓ and
C^n,m_ℓ, respectively, as ℓ→∞,
and these limits form a solution pair for the Galerkin approximation (<ref>), (<ref>).
For further details, see <cit.>.  This establishes the existence of a solution to the Galerkin approximations
(<ref>), (<ref>) for any fixed pair of integers n,m∈ℕ.
The existence of a discrete solution triple for (<ref>)–(<ref>) then follows by the discrete inf-sup condition
stated in Proposition <ref>, and we write P^n,m=∑^Ñ_n_i=1γ^n,m_iy_i
where {y_i}^Ñ_n_i=1 is a basis of ^n_0.

We are now ready to state and prove our main theorem in this section. It asserts that, as n,m→∞, the sequence of discrete
solution triples converges to a weak solution triple of the regularized problem.


Suppose that Ω⊂^3 is a convex polyhedral domain and c_d∈ W^1,s(Ω) for some s>3.
Let us assume that r:_≥0→_≥0 is a Hölder-continuous function with
r^->3/2>t/t-2, t>2, and let f∈ (W^1,r^-_0(Ω)^3)^*.
Let (^n,m,P^n,m,C^n,m)∈^n_div×^n_0×(^m+c_d) be a discrete solution
triple defined by the finite element approximation (<ref>)–(<ref>).
Then, the following convergence results hold.



  * At the first level of Galerkin approximation, there exists a subsequence (not relabelled) with respect to m such that (as m→∞),

    ^n,m   →^n      ,
    ^n,m   →^n      ,
    
    P^n,m   → P^n      ,
    
    C^n,m   ⇀ C^n      weakly  in  W^1,2(Ω),

where ^n∈^n, P^n∈^n_0.

  * At the second level of Galerkin approximation, there exists a subsequence (not relabelled) with respect to n such that (as n→∞),

    ^n    ⇀      weakly  in  W^1,r^-_0(Ω)^3
    
    P^n    ⇀ p       weakly  in  L^j'(Ω)    ∀  j>max{r^+,2},
    
    C^n    ⇀ c       weakly  in  W^1,2(Ω),
    
    C^n    → c       strongly  in  C^0,α(Ω),

where (,p,c)=(^k,p^k,c^k) is a weak solution triple of the regularized problem (<ref>)–(<ref>).





Proof of Theorem <ref>
The limit m→∞
First, we shall derive some uniform bounds, independent of m∈ℕ, and let m tend
to infinity by using the weak compactness properties in the corresponding reflexive spaces. For simplicity, we shall denote ^n,m(C^n,m,^n,m), _c^n,m_c(C^n,m,∇ C^n,m,^n,m).

We test with ^n,m∈^n_div in (<ref>); then, thanks to the skew symmetry of B_u[·,·,·], we have


    ∫_Ω^n,m·∇^n,m+1/k|^n,m|^t=∫_Ω^n,m·^n,m+1/k|^n,m|^t=⟨f,^n,m⟩.


By (<ref>) and Young's inequality, we have

    ∫_Ω|∇^n,m|^r(C^n,m)+|^n,m|^r'(C^n,m)+|^n,m|^t≤ C_1,

where C_1 is independent of m.


Next, we test with C^n,m-c_d ∈^m in (<ref>) and deduce that


    ∫_Ω_c(C^n,m,∇ C^n,m^n,m)·∇(C^n,m-c_d)=B_c[C^n,m,^n,m,c_d].


By (<ref>), (<ref>), Hölder's inequality and Young's inequality,

    ∇ C^n,m^2_2
       ≤∫_Ω|∇ C^n,m||∇ c_d|+B_c[C^n,m,^n,m,c_d]
       ≤ε∇ C^n,m^2_2+C(ε)∇ c_d^2_2+B_c[C^n,m,^n,m,c_d].

Then, by Sobolev embedding,

    B_c[C^n,m,^n,m,c_d]
       =1/2∫_Ωc_d^n,m·∇ C^n,m-1/2∫_ΩC^n,m^n,m·∇ c_d
       =∫_Ωc_d^n,m·∇ C^n,m+1/2∫_ΩC^n,m(div ^n,m)c_d
       ≤c_d_∞^n,m_2∇ C^n,m_2+1/2 c_d_∞ C^n,m_(r^-)'div ^n,m_r^-
       ≤ C^n,m_1,r^-∇ C^n,m_2+C^n,m_1,r^-∇ C^n,m_3r^-/4r^–3
       ≤ C(ε)^n,m^2_1,r^-+ε∇ C^n,m^2_2.

Hence, by (<ref>) and (<ref>), we have

    ∫_Ω|∇ C^n,m|^2+|^n,m_c|^2≤ C_2,

where C_2 is independent of m.

Next, we shall derive a uniform bound on the pressure. By Proposition <ref> together with (<ref>), (<ref>)
and the equivalence of norms in the finite-dimensional spaces, we have

    β_rP^n,m_(r^+)'   ≤sup_0≠∈^n⟨div ,P^n,m⟩/_1,r^+
       ≤sup_0≠∈^n|∫_Ω^n,m·|/_1,r^++Csup_0≠∈^n|B_u[^n,m,^n,m,]-⟨f,⟩|/_1,r^-
       ≤ C sup_0≠∈^n^n,m_(r^+)'_r^+/_1,r^++C(n)sup_0≠∈^n^n,m^2_2(r^-)'_1,r^-+f_-1_1,r^-/_1,r^-.

Therefore, by (<ref>), we deduce that

    P^n,m_(r^+)'≤ C(n).


Now we are ready to let m tend to infinity. By (<ref>) and (<ref>) with the equivalence of norms in
finite-dimensional spaces, we have |α^n,m|≤ C(n) and |γ^n,m|≤ C(n).
Then, together with the uniform estimates (<ref>), we can extract (not relabelled) subsequences such that

    α^n,m   →α^n      strongly   in   ^N_n,
    γ^n,m   →γ^n      strongly   in   ^Ñ_n,
    
    C^n,m   ⇀ C^n      weakly   in   W^1,2(Ω).

From (<ref>), (<ref>) and compact embedding,  we have

    ^n,m   →^n      ,
    ^n,m   →^n      ,
    
    P^n,m   → P^n      ,
    
    C^n,m   → C^n      strongly   in   L^2(Ω).

By (<ref>) and (<ref>), note that

    ^n∈^n    and     P^n∈^n_0.

Finally, from (<ref>), we can extract a further subsequence (not relabelled) such that

    C^n,m→ C^n    a.e.   in   Ω.


Note that since  is continuous, by (<ref>) and (<ref>), we have

    (C^n,m,^n,m)→(C^n,^n)    a.e.   in   Ω.

Now, by (<ref>), we have that, for sufficiently large m∈ℕ,

    |^n,m|<1+|^n|    .

Thus, by (<ref>), we have, for sufficiently large m∈ℕ,

    |(C^n,m,^n,m)|
       ≤ C|^n,m|^r(C^n,m)-1+C
       ≤ C(1+|^n|)^r(c^n,m)-1+C
       ≤ C(1+|^n|)^r^+-1+C,

and C(1+|^n|)^r^+-1+C∈ L^(r^+)'(Ω). Therefore, by the Dominated Convergence Theorem, we have

    ^n,m→^n(C^n,^n)    strongly   in   L^(r^+)'(Ω)^3×3.


Furthermore, by (<ref>) and (<ref>), together with the Dominated Convergence Theorem,

    K(C^n,m,|^n,m|)→K(C^n,|^n|)    strongly   in   L^q(Ω)    ∀  q∈(1,∞).

Therefore, together with (<ref>), we have

    ^n,m_c⇀^n_c_c(C^n,∇ C^n,^n)    weakly   in   L^2(Ω)^3.


Now we are ready to pass m to infinity in the Galerkin approximation (<ref>)–(<ref>). First, by (<ref>) and (<ref>),

    B_u[^n,m,^n,m,]   → B_u[^n,^n,]      ∀ ∈^n,
    1/k|^n,m|^t-2^n,m·   →1/k|^n|^t-2^n·      ∀ ∈^n.


Furthermore, from (<ref>) and (<ref>),

    ∫_Ω^n,m·   →∫_Ω^n·      ∀ ∈^n,
    ⟨div ,P^n,m⟩   →⟨div ,P^n⟩      ∀ ∈^n.

Therefore, we have

    ∫_Ω^n·+1/k|^n|^t-2^n·+B_u[^n,^n,]-⟨div ,P^n⟩=⟨f,⟩    ∀ ∈^n.

Moreover, from (<ref>) and (<ref>),

    ∫_ΩQ div ^n=0    ∀  Q∈^n.



Next, let us investigate the limit of the concentration equation, (<ref>). We fix an arbitrary Z∈ W^1,2_0(Ω)
and define Z^mΠ^m_Z∈^m. Thanks to (<ref>) and (<ref>),

    C^n,m^n,m-C^n^n_2
       ≤(^n,m-^n)C^n,m_2+^n(C^n,m-C^n)_2
       ≤^n,m-^n_∞C^n,m_2+^n_∞C^n,m-C^n_2→0.

Also, thanks to (<ref>) and (<ref>),

    Z^m^n,m-Z^n_2
       ≤(^n,m-^n)Z^m_2+^n(Z^m-Z)_2
       ≤^n,m-^n_∞Z^m_2+^n_∞Z^m-Z_2→0.

In other words, we have

    C^n,m^n,m   → C^n^n      strongly   in   L^2(Ω)^3,
    
    Z^m^n,m   → Z^n      strongly   in   L^2(Ω)^3.

By (<ref>) and (<ref>),

    |∫_ΩZ^m^n,m·∇ C^n,m-∫_ΩZ^n·∇ C^n|
         ≤∫_Ω|Z^m^n,m-Z^n||∇ C^n,m|+|∫_ΩZ^n(∇ C^n,m-∇ C^n)|→0.

Moreover, from (<ref>) and (<ref>),

    |∫_ΩC^n,m^n,m∇ Z^m-∫_ΩC^n^n·∇ Z|
         ≤C^n,m^n,m_2Z^m-Z_1,2+Z_1,2C^n,m^n,m-C^n^n_2→0.

Therefore, we have

    lim_m→∞B_c[C^n,m,^n,m,Z^m]=B_c[C^n,^n,Z].

Finally, from (<ref>),

    ∫_Ω^n,m_c·∇ Z^m→∫_Ω^n_c·∇ Z    as   m→∞.

Altogether, we have

    ∫_Ω^n_c·∇ Z+B_c[C^n,^n,Z]=0    ∀  Z∈ W^1,2_0(Ω).



The limit n→∞
Now we shall derive further uniform estimates and let n pass to infinity. First, we test with ^n in (<ref>).
Then, by (<ref>) and (<ref>), we have

    ∫_Ω^n·^n+1/k|^n|^t=⟨f,^n⟩.

By using (<ref>) and Young's inequality, we have

    ∫_Ω|^n|^r(C^n)+|^n|^r'(C^n)+1/k|^n|^t≤ C_1,

where C_1 is independent of n, which leads us to

    ^n_1,r^-^r^-+^n_(r^+)'^(r^+)'+1/k^n^t_t≤ C_1,

where C_1 is independent of n.

Next, we test with C^n-c_d in (<ref>), and by (<ref>) we obtain

    ∫_Ω^n_c·∇ C^n=∫_Ω^n_c·∇ c_d+B_c[C^n,^n,c_d].

From (<ref>), (<ref>), Hölder's inequality and Young's inequality we have

    ∇ C^n^2_2
       ≤ C∫_Ω|∇ C^n||∇ c_d|+B_c[C^n,^n,c_d]
       ≤ε∇ C^n^2_2+C(ε)∇ c_d^2_2+B_c[C^n,^n,c_d].

Furthermore, by Sobolev embedding,

    B_c[C^n,^n,c_d]
       =1/2∫_Ωc_d^n·∇ C^n-1/2∫_ΩC^n^n·∇ c_d
       =∫_Ωc_d^n·∇ C^n+1/2∫_ΩC^n(div ^n)c_d
       ≤c_d_∞^n_2∇ C^n_2+c_d_∞/2C^n_(r^-)'div ^n_r^-
       ≤ C^n_1,r^-∇ C^n_2+C^n_1,r^-∇ C^n_3r^-/4r^–3
       ≤ C(ε)^n^2_1,r^-+ε∇ C^n^2_2.

Hence, from (<ref>) and (<ref>),

    ∫_Ω|∇ C^n|^2+|_c^n|^2≤ C_2,

where C_2 is independent of n. Thus we have

    C^n^2_1,2+^n_c^2_2≤ C_2,

where C_2 is independent of n.

Now, since 3/2>t/t-2, by Sobolev embedding and the uniform estimates (<ref>) and (<ref>), for s>3
sufficiently close to 3,


    C^n^n_s≤C^n_6^n_6s/6-s≤ CC^n_1,2^n_t≤ C,


where C is independent of n.

Also, for s>3 sufficiently close to 3, we have

    ∇ C^n·^n_3s/s+3≤∇ C^n_2^n_6s/6-s≤ CC^n_1,2^n_t≤ C,


where C is independent of n.

Therefore, we can apply Theorem <ref> with F=C^n^n and g=∇ C^n·^n. Hence, there exists
an α_1∈(0,1) such that

    C^n_C^0,α_1(Ω)≤ C_3.

Since C^0,α_1(Ω)↪↪ C^0,α̃_1(Ω) for all α̃_1∈(0,α_1), we have

    C^n→ c    strongly   in   C^0,α̃_1(Ω),

which implies that

    r∘ C^n→ r∘ c    strongly   in   C^0,β_1(Ω),

for some β_1∈(0,1).

We now apply Proposition <ref>. For a given r^+>0, choose j>max{r^+,2}. Then, since r^->3/2,
we have that W^1,j_0(Ω)^3↪ L^2(r^-)'(Ω)^3 by Sobolev embedding.
Furthermore, since t/t-2<r^-, we have that 2(r^-)'<t. Now, from (<ref>) and (<ref>),

    β_rP^n_j'   ≤sup_0≠∈^n⟨div ,P^n⟩/_1,j
       ≤sup_0≠∈^n|∫_Ω^n·+B_u[^n,^n,]-⟨f,⟩|/_1,j
       ≤ Csup_0≠∈^n^n_(r^+)'_1,r^+/_1,r^+
       +Csup_0≠∈^n^n^2_t_1,r^-+f_-1_1,r^-/_1,r^-
       +Csup_0≠∈^n^n_2(r-)'_2(r-)'^n_1,r^-/_2(r-)'.

Hence, by noting (<ref>),

    P^n_j'≤ C_4,

where C_4 is independent of n.

Now, by (<ref>)–(<ref>), thanks to the reflexivity of the relevant spaces
and by compact embedding, we can extract (not relabelled) subsequences such that


    ^n    ⇀      weakly  in  W^1,r^-_0(Ω)^3∩ L^t(Ω)^3, 
    ^n    →      strongly  in  L^σ(Ω)^3    ∀ σ∈[1,t),
    
    |^n|^t-2^n   ⇀||^t-2      weakly  in  L^t/t-1(Ω)^3,
    
    C^n    ⇀ c       weakly  in  W^1,2(Ω),
    
    C^n    → c       strongly  in  C^0,α̃_1(Ω),
    
    P^n    ⇀ p       weakly  in  L^j'(Ω)    ∀  j>max{r^+,2},
    ^n    ⇀      weakly  in  L^(r^+)'(Ω)^3× 3,
    _c^n    ⇀_c       weakly  in  L^2(Ω)^3.

Before proceeding further, we note that these limits, together with weak lower semicontinuity and  (<ref>),
in conjunction with Korn's inequality, imply that

    ∫_Ω|∇|^r(c)+||^r'(c)≤ C,

hence the limit function  is, in fact, contained in the space W^1,r(c)_0(Ω)^3; see <cit.> for the
details of the proof of this.

Next, we shall prove that the limit function  is pointwise divergence-free. For an arbitrary q∈ C^∞_0(Ω),
by (<ref>),

    0   =∫_Ω(Π^n_q) div ^n
       =∫_Ω(Π^n_q-q) div ^n+∫_Ωq(div ^n-div )+∫_Ωq div .

The first term tends to zero by (<ref>), (<ref>) and the second term converges to zero by (<ref>). Therefore,

    ∫_Ωq div =0    for   any   q∈ C^∞_0(Ω),

which implies that div  = 0 a.e. on Ω.

Now, we shall identify the limit of the convective term B_u[·,·,·] as follows. For an arbitrary ∈ W^1,∞_0(Ω)^3, we define ^nΠ^n_div∈^n. Then, by (<ref>), we have

    ^n→    strongly   in   W^1,σ_0(Ω)^2   for   σ∈(1,∞).

By (<ref>),

    ^n⊗^n→⊗    strongly   in   L^1+ε(Ω)^3×3.

Hence, we can identify the second part of the convective term

    -∫_Ω(^n⊗^n)·∇^n→-∫_Ω(⊗)·∇    as   n→∞.


Also, we assert that ^n·^n→· strongly in L^(r^-)'(Ω). Indeed,

    ^n·^n-·_(r^-)'   ≤(^n-)^n+(^n-)_(r^-)'
       ≤^n-_σ^n_t-ε+^n-_t-ε_σ

for some σ∈(1,∞). The first term tends to zero thanks to (<ref>), (<ref>) and the second term tends to zero by (<ref>). Therefore, since div =0, we have

    ∫_Ω(^n⊗^n)·∇^n   =-∫_Ω(^n⊗^n)·∇^n+∫_Ω(div ^n) ^n·^n
       →-∫_Ω(⊗)·∇    as   n→∞.

Altogether, we then deduce that

    lim_n→∞B_u[^n,^n,^n]=-∫_Ω(⊗)·∇.


Now, we are ready to pass n to infinity in the Navier–Stokes equations. Since Π^n_div is linear, by noting (<ref>), we have

    ⟨div ,P^n⟩   =⟨div ^n,P^n⟩+⟨div (-^n),P^n⟩
       =∫_Ω(C^n,^n)·^n+1/k|^n|^t-2^n·^n -⟨f,^n⟩+B_u[^n,^n,^n]
             +⟨div (-^n),P^n⟩
       →∫_Ω·+1/k||^t-2·+div(⊗)·-⟨f,⟩,

where we have used (<ref>), (<ref>), (<ref>), (<ref>) and (<ref>). Also, by (<ref>) again,

    ⟨div ,P^n⟩→⟨div ,p⟩.

Collecting all the limits gives us

    ∫_Ω·+1/k||^t-2·+div (⊗)·-⟨div ,p⟩=⟨f,⟩    ∀ ∈ W^1,∞_0(Ω)^3.


With the same argument as above, we also have that

    ∫_Ω·+1/k||^t-2·+div (⊗)·=⟨f,⟩    ∀ ∈ W^1,∞_0,div(Ω)^3.


Note that by Proposition <ref> and (<ref>), we have

    p∈ L^r'(c)_0(Ω).


Now, let us investigate the limit of the convection-diffusion equation, (<ref>). For an arbitrary but fixed z∈ W^1,2_0(Ω),
we define Z^nΠ^n_z∈^n. Thanks to (<ref>) and (<ref>),

    C^n^n-c_2
       ≤(C^n-c)^n_2+c(^n-)_2
       ≤C^n-c_∞^n_2+c_∞^n-_2→0.

Moreover, by (<ref>), (<ref>) and Sobolev embedding,

    Z^n^n-z_2
       ≤(Z^n-z)^n_2+z(^n-)_2
       ≤Z^n-z_6^n_3
    +z_6^n-_3
       ≤ CZ^n-z_1,2^n_3+Cz_1,2^n-_3→0.

In other words,

    C^n^n   → c    strongly   in   L^2(Ω)^3,
    
    Z^n^n   → z    strongly   in   L^2(Ω)^3.

From (<ref>) and (<ref>),

    |∫_ΩZ^n^n·∇ C^n-∫_Ωz·∇ c|
           ≤∫_Ω|Z^n^n-z||∇ C^n|+|∫_Ωz·(∇ C^n-∇ c)|→0.

Therefore, as div =0 a.e. on Ω, we obtain

    ∫_ΩZ^n^n·∇ C^n→∫_Ωz·∇ c=-∫_Ωc·∇ z    as   n→∞.

Moreover, by (<ref>) and (<ref>),

    |∫_ΩC^n^n·∇ Z^n-∫_Ωc·∇ z|
           ≤C^n^n_2Z^n-z_1,2+C^n^n-c_2z_1,2→0.

Altogether, we have

    lim_n→∞ B_c[C^n,^n,Z^n]=-∫_Ωc·∇ z.

Finally, by (<ref>) and (<ref>), we have

    ∫_Ω_c(C^n,∇ C^n,^n)·∇ Z^n→∫_Ω_c·∇ z    as   n→∞.

By collecting all the limits, we obtain that

    ∫_Ω_c·∇ z-c·∇ z=0    ∀  z∈ W^1,2_0(Ω).

As we can see from (<ref>) and (<ref>), what we now need to prove is the identification of the limits:

    =(c,)     and     _c=_c(c,∇ c, ).

To this end, we require the following lemma.



The sequences {^n}_n ∈ℕ and {C^n}_n ∈ℕ satisfy the following equality:

    lim_n→∞∫_Ω(((C^n,^n)-(C^n,))·(^n-))^1/4=0.



The detailed proof of Lemma <ref> is presented in Section 4.2 in <cit.>. Here, we shall briefly summarize the key steps of the proof as we shall require a similar, but more involved, argument in the next section. The strategy is to decompose the integral into several terms and to estimate them separately. To this end, for arbitrary but fixed χ>0, we introduce the matrix-truncation function T_χ:^3× 3→^3× 3 by


T_χ(M)={ [          M for |M|≤χ,;     χM/|M| for |M|>χ. ].


The essential step in the proof relies on using a discrete Lipschitz truncation technique. In  <cit.> a version of the discrete Lipschitz truncation method in variable-exponent norms was presented: see Theorem 3.15 in <cit.> and let ^n_j denote the discrete Lipschitz truncation of the function of ^n.

The most important and difficult part of the proof is to estimate the following term:

    lim_χ→ 0lim_j →∞lim_n →∞∫_Ω((C^n,^n)-(C^n,T_χ()))·(^n_j-T_χ())≤ 0;

(see eq. (4.23) in <cit.>).
The other terms arising from the decomposition can be easily estimated by using the uniform bound (<ref>), Hölder's inequality and the discrete Lipschitz truncation theorem, Theorem 3.15 in <cit.>.

To estimate (<ref>), we introduce the following discretely divergence-free approximations with zero trace on ∂Ω:

    Ψ^n_j   ℬ^n(div ^n_j),
    Φ^n_j   ^n_j-Ψ^n_j.

Here ℬ^n is a discrete Bogovskiĭ operator defined in Section 3.4 of <cit.>. It is then clear that Φ^n_j has zero trace on ∂Ω and, by construction, Φ^n_j∈^n_div. Moreover, it can be easily verified, by using basic properties of the discrete Lipschitz truncation and the discrete Bogovskiĭ operator, that

    Φ^n_j   ⇀_j-ℬ(div _j)Φ_j      weakly  in   W^1,σ_0(Ω)^3,
    Φ^n_j   →Φ_j      strongly  in  L^σ(Ω)^3,

as n→∞, where σ∈(1,∞) is arbitrary. We can then rewrite (<ref>) above in terms of this approximation to obtain

    ∫_Ω((C^n,^n)-(C^n,T_χ()))·(^n_j-T_χ())
       =∫_Ω(C^n,^n)·(Φ^n_j+Ψ^n_j)
             -∫_Ω(C^n,^n)· T_χ()-∫_Ω(C^n,T_χ())·(^n_j-T_χ())
        B^n,1_χ,j-B^n,2_χ,j-B^n,3_χ,j.

Now we use (<ref>) with =Φ^n_j∈^n_div and pass to the limit; thus we have, by (<ref>), that

    lim_n→∞∫_Ω^n·Φ^n_j   =-lim_n→∞B_u[^n,^n,Φ^n_j]-∫_Ω1/k|^n|^t-2^n·Φ^n_j+lim_n→∞⟨f,Φ^n_j⟩
       =∫_Ω(⊗)·∇Φ_j-1/k||^t-2·Φ_j+⟨f,Φ_j⟩
       =∫_Ω·Φ_j.


Furthermore, with the help of Lipschitz truncation, we can show that

    lim_n→∞∫_Ω^n·Ψ^n_j   ≤(C/2^j/r^+)^γ(r^-,r^+),
    ∫_Ω·ℬ(div _j)   ≤(C/2^j/r^+)^γ(r^-,r^+).

Altogether, we have

    lim_χ→∞lim_j→∞lim_n→∞(B^n,1_χ,j-B^n,2_χ,j-B^n,3_χ,j)≤lim_χ→∞∫_Ω(-(c,T_χ()))·(_j-T_χ())

The last limit is equal to zero by using the Dominated Convergence Theorem. That completes the proof of (<ref>), and thereby also of the most technical step in the proof of the lemma.



 Now we are ready to identify the limits. In the above lemma, since the integrand is nonnegative, (<ref>) also holds with Ω replaced by the set
Q_γ⊂Ω defined by

    Q_γ{x∈Ω:||≤γ},

with a given  γ>0; thus, from the sequence of integrands featuring in  (<ref>), we can extract a subsequence
(again not relabelled), which converges to zero almost everywhere in Q_γ. Then, by Egoroff's Theorem,
for an arbitrary ε>0, there exists a subset Q^ε_γ⊂ Q_γ⊂Ω satisfying
|Q_γ∖ Q^ε_γ|<ε, where the convergence of  integrands is uniform.
Note that, thanks to the choice of Q^ε_γ, we have


    lim_γ→∞lim_ε→0|Ω∖ Q^ε_γ|=
    lim_γ→∞lim_ε→0[ |Ω∖ Q_γ| + |Q_γ∖ Q^ε_γ|
    ]= 0.


Moreover, we have from the uniform convergence of the integrands that

    lim_n→∞∫_Q^ε_γ((C^n,^n)-(C^n,))· (^n-)=0.

Since  is bounded on Q^ε_γ, by the Dominated Convergence Theorem we have (C^n,)→(c,) strongly in L^q(Ω)^3× 3 for any q∈[1,∞).
Hence, from the above L^q-convergence, (<ref>), and (<ref>), we obtain

    lim_n→∞∫_Q^ε_γ(C^n,^n)·(^n-)=0.

Thus, by the boundedness of  on Q^ε_γ and  (<ref>), we have

    lim_n→∞∫_Q^ε_γ(C^n,^n)·^n=∫_Q^ε_γ·.

Now, let B∈ L^∞(Q^ε_γ)^3× 3 be arbitrary but fixed.
From the monotonicity (<ref>), (<ref>), the L^q-convergence of (C^n,B)→(c,B) and the weak convergence (<ref>), we have

    0   ≤lim_n→∞∫_Q^ε_γ((C^n,^n)-(C^n,B))·(^n-B)
       =∫_Q^ε_γ·(-B)-∫_Q^ε_γ(c,B)·(-B)
       =∫_Q^ε_γ(-(c,B))·(-B).

Now we are ready to use Minty's trick. First, we choose B=±λA with λ>0 and A∈ L^∞(Q^ε_γ)^3× 3. Then, passing to the limit λ→0,
the continuity of  gives us


    ∫_Q^ε_γ(-(c,))·A=0.

Hence, we have that

    =(c,)   a.e.   on   Q^ε_γ.

Now we pass ε→0 and then γ→∞ to conclude that

    =(c,)   a.e.   on   Ω.


Finally, since  is strictly monotonic and C^n→ c in C^0,α̃_1(Ω),
by (<ref>) we deduce that

    ^n→   a.e.   on   Ω.

By the Dominated Convergence Theorem, with (<ref>), (<ref>) and (<ref>), we obtain that

    _c(C^n,∇ C^n,^n)⇀_c(c,∇ c, )    weakly  in  L^2(Ω)^3.

Therefore, by the uniqueness of the weak limit, we can identify

    _c=_c(c,∇ c,).





Proof of Theorem <ref>
Minimum and maximum principles
Before we proceed, let us prove minimum and maximum principles for the concentration. Let φ^k_1=(c^k-min_x∈∂Ωc_d)_-
and φ^k_2=(c^k-max_x∈∂Ωc_d)_+. Since c^k=c_d on ∂Ω, it is clear that
φ^k_1,φ^k_2∈ W^1,2_0(Ω), so we can test with φ^k_1 and φ^k_2 in (<ref>). Therefore, we have

    -∫_Ω^kc^k·∇φ^k_1+∫_Ω_c∇φ^k_1=0,


    -∫_Ω^kc^k·∇φ^k_2+∫_Ω_c∇φ^k_2=0.

We first consider (<ref>). From (<ref>) with integration by parts we obtain

    ∫_Ω^-^k·∇ c^kφ^k_1+∫_Ω^-C|∇ c^k|^2≤0,

where Ω^-={x∈Ω:φ^k_1(x)<0}, since div ^k=0 and ^k=0 on ∂Ω. By using the fact that ∇ c^k=∇φ^k_1 on Ω^- and the extension of ∇ c^k from Ω^- to the whole domain Ω by using the negative part, we have

    ∫_Ω^k·∇φ^k_1 φ^k_1+∫_ΩC|∇φ^k_1|^2≤0.

Note that

    ∫_Ω^k·∇φ^k_1φ^k_1=1/2∫_Ω^k·∇|φ^k_1|^2=-1/2∫_Ω(div ^k)|φ^k_1|^2=0,

and thus,

    φ^k_1=(c^k-min_x∈∂Ωc_d)_-=constant  a.e.  in  Ω.

In the same way, we can also show that

    φ^k_2=(c^k-max_x∈∂Ωc_d)_+=constant  a.e.  in  Ω.

By combining the above results we finally obtain that

    min_x∈∂Ωc_d≤ c^k≤max_x∈∂Ωc_d    a.e.   in   Ω.



The limit k→∞
First, note that by weak lower semicontinuity of the norm-function, and (<ref>), (<ref>) and (<ref>), we obtain the following
uniform estimates, independent of k∈ℕ:

    ^k^r^-_1,r^-+(c^k,^k)^(r^+)'_(r^+)'+1/k^k^t_t≤ C_1,


    c^k^2_1,2+_c(c^k,∇ c^k,^k)^2_2≤ C_2,


    p^k^j'_j'≤ C_3,

for some positive constants C_1, C_2 and C_3, which are independent of k∈ℕ.

Now, since r^->3/2, by the min/max principle (<ref>), Sobolev embedding and the uniform estimate (<ref>),
for s>3 sufficiently close to 3,

    c^k^k_s≤c^k_∞^k_s≤ C^k_1,r^-≤ C.


Therefore, we can again apply Theorem <ref> with F=c^k^k and g=0. Hence, there exists an α_2∈(0,1) such that

    c^k_C^0,α_2(Ω)≤ C_4,

for some positive constant C_4 independent of k∈ℕ. Since C^0,α_2(Ω)↪↪ C^0,α̃_2(Ω) for all α̃_2∈(0,α_2), we have

    c^k→ c    strongly   in   C^0,α̃_2(Ω),

which implies that

    r∘ c^k→ r∘ c    strongly   in   C^0,β_2(Ω),

for some β_2∈(0,1).

Therefore, by the reflexivity of the relevant spaces and compact embedding, there exists a subsequence (not relabelled) such that

    ^k    ⇀      weakly  in  W^1,r^-_0,div(Ω)^3, 
    ^k    →      strongly  in  L^2(1+ε)(Ω)^3,
    
    c^k    ⇀ c       weakly  in  W^1,2(Ω),
    
    c^k    → c       strongly  in  C^0,α̃_2(Ω),
    
    p^k    ⇀ p       weakly  in  L^j'(Ω)    ∀  j>max{r^+,2},
    (c^k,^k)    ⇀      weakly  in  L^(r^+)'(Ω)^3× 3,
    _c(c^k,∇ c^k, ^k)    ⇀_c       weakly  in  L^2(Ω)^3.


Again, by the weak lower semicontinuity of norms, (<ref>) and (<ref>) together with Korn's inequality, we have that

    ∫_Ω|∇|^r(c)+||^r'(c)≤ C,

and thus the weak solution  is in the desired space W^1,r(c)_0(Ω)^3.

Now we shall let k→∞ in (<ref>), with ∈ W^1,∞_0(Ω)^3 chosen arbitrarily. By (<ref>),

    ^k⊗^k→⊗    strongly   in   L^1+ε(Ω)^3×3.

Thus, we can identify the limit of the convective term

    -∫_Ω(^k⊗^k)·∇→-∫_Ω(⊗)·∇    as   k→∞,    ∀ ∈ W^1,∞_0(Ω)^3.


Next, by (<ref>), we have that

    1/k^k^t-1_t→0    as   k→∞.

Therefore, we have

    1/k|∫_Ω|^k|^t-2^k·|≤1/k^k^t-1_t_t
    →0    as   k→∞,    ∀ ∈ W^1,∞_0(Ω)^3.


We recall from the identification asserted in (<ref>) that =(c,) a.e. on Ω; more precisely, with the index k reinstated in our notation, ^k=(c^k,^k) a.e. on Ω. Hence, from (<ref>) and (<ref>), we obtain

    ⟨div ,p^k⟩→⟨div ,p⟩  and  ∫_Ω^k·→∫_Ω·  as   k→∞,  ∀ ∈ W^1,∞_0(Ω)^3.


Altogether, we have

    ∫_Ω·+(⊗)·∇-⟨div ,p⟩=⟨f,⟩    ∀ ∈ W^1,∞_0(Ω)^3.


Furthermore, it is clear that

    ∫_Ω·+(⊗)·∇=⟨f,⟩    ∀ ∈ W^1,∞_0,div(Ω)^3.


Note that by Proposition <ref> and (<ref>) we have

    p∈ L^r'(c)_0(Ω).


Now, let us investigate the limit of the concentration equation (<ref>).
Let us choose an arbitrary, but fixed, z∈ W^1,2_0(Ω). By (<ref>) and (<ref>),

    c^k^k-c_2
    ≤(c^k-c)^k_2+c(^k-)_2
    ≤c^k-c_∞^k_2+c_∞^k-_2→0.

In other words,

    c^k^k→ c    strongly   in   L^2(Ω)^3.

Hence we have

    ∫_Ωc^k^k·∇ z→∫_Ωc·∇ z.


Recalling the identification (<ref>) and reinstating the index k, we have
_c^k:=_c(c^k,∇ c^k, ^k); hence, by (<ref>), we get


    ∫_Ω_c^k·∇ z→∫_Ω_c·∇ z    as   k→∞.


By collecting the above limits, we deduce that

    ∫_Ω_c·∇ z-c·∇ z=0    ∀  z∈ W^1,2_0(Ω).


As a final step, we need to identify the limits:

    =(c,)     and     _c=_c(c,∇ c, ).

To this end, analogously as before, we need to prove the following equality:

    lim_k→∞∫_Ω(((c^k,^k)-(c^k,))·(^k-))^1/4=0.

The proof is similar to the one presented in the previous section. The only part of the argument that we shall give here in detail is the proof of the analogue of (<ref>)–(<ref>) since we now have a different weak formulation at this level. The other parts of the proof proceed as Section 4.2 in <cit.>.

First we define a divergence-free approximation with zero trace as follows:

    Φ^k_j^k_j-ℬ(div ^k_j),

where ℬ is the Bogovskiĭ operator introduced in Theorem <ref>. Then, as before, we have

    Φ^k_j   ⇀_j-ℬ(div _j)Φ_j      weakly  in   W^1,σ_0(Ω)^3,
    Φ^k_j   →Φ_j      strongly  in  L^σ(Ω)^3,

as k→∞, where σ∈(1,∞) is arbitrary.

Let us further define χ^n,k_1,jΠ^n_divΦ^k_j. Then, by (<ref>),

    χ^n,k_1,j→Φ^k_j    strongly   in   W^1,σ_0(Ω)^3,   ∀ σ∈(1,∞).

Now, by (<ref>),

    ∫_Ω^n·χ^n,k_1,j=-B_u[^n,^n,χ^n,k_1,j]-∫_Ω1/k|^n|^t-2^n·χ^n,k_1,j+⟨f,χ^n,k_1,j⟩.

If we take n→∞ in the above equality, we have

    ∫_Ω(c^k,^k)·Φ^k_j=∫_Ω(^k⊗^k)·∇Φ^k_j-1/k|^k|^t-2^k·Φ^k_j+⟨f,Φ^k_j⟩.


Next, we define χ^n,k_2,jΠ^n_divΦ_j, and then we have

    χ^n,k_2,j→Φ_j    strongly   in   W^1,σ_0(Ω)^3,   ∀ σ∈(1,∞).

Again, by (<ref>),

    ∫_Ω^n·χ^n,k_2,j=-B_u[^n,^n,χ^n,k_2,j]-∫_Ω1/k|^n|^t-2^n·χ^n,k_2,j+⟨f,χ^n,k_2,j⟩.

If we take n→∞, we have

    ∫_Ω(c^k,^k)·Φ_j=∫_Ω(^k⊗^k)·∇Φ_j-1/k|^k|^t-2^k·Φ_j+⟨f,Φ_j⟩.

Subsequently, if we pass k to the infinity, we obtain

    ∫_Ω·=∫_Ω(⊗)·∇Φ_j+⟨f,Φ_j⟩.


Therefore, from (<ref>) and (<ref>), we deduce that

    lim_k→∞∫_Ω(c^k,^k)·Φ^k_j   =lim_k→∞∫_Ω(^k⊗^k)·∇Φ^k_j-1/k|^k|^t-2^k·Φ^k_j+lim_k→∞⟨f,Φ^k_j⟩
       =∫_Ω(⊗)·∇Φ_j+⟨f,Φ_j⟩
       =∫_Ω·Φ_j,

which is the desired analogue of (<ref>)–(<ref>) corresponding to the limit k→∞, and thereby the proof of (<ref>) has been completed.

We can then use the same argument as the one we employed in the previous section to
identify = ^k=(c^k,^k) and _c=_c^k = _c(c^k,∇ c^k,^k)
(cf. (<ref>) and (<ref>), with the index k reinstated),
and thus we can again identify =(c,), _c=_c(c,∇ c,).
That completes the proof of the convergence theorem.





§ CONCLUSIONS


We have considered a system of nonlinear partial differential equations modelling the motion of an incompressible chemically reacting generalized Newtonian fluid in three space dimensions. The governing system consists of a steady convection-diffusion equation for the concentration and a generalized steady power-law-type fluid flow model for the velocity and the pressure,  where the viscosity depends on both the shear-rate and the concentration through a concentration-dependent power-law index.
We performed a rigorous convergence analysis of a finite element approximation of a regularized counterpart of the model;
specifically, we showed the convergence of the finite element method to a weak solution of the regularized model. We then
proved that weak solutions of the regularized problem converge to a weak solution of the original problem.



§ ACKNOWLEDGEMENTS

Seungchan Ko's work was supported by the UK Engineering and Physical Sciences Research Council [EP/L015811/1].



abbrv





