A virtual element method for a convective Brinkman-Forchheimer problem coupled with a heat equation
Danilo Amigo Felipe Lepe Enrique Otarola Gonzalo Rivera


1/f Noise in the Heliosphere: A Target for PUNCH Science
    [
    September 9, 2024
========================================================



§ ABSTRACT 
We develop a virtual element method to solve a convective Brinkman-Forchheimer problem coupled with a heat equation. This coupled model may allow for thermal diffusion and viscosity as a function of temperature. Under standard discretization assumptions, we prove the well posedness of the proposed numerical scheme. We also derive optimal error estimates under appropriate regularity assumptions for the solution. We conclude with a series of numerical tests performed with different mesh families that complement our theoretical findings.




non-isothermal flows, nonlinear equations, a convective Brinkman-Forchheimer problem, a heat equation, virtual element methods, stability, a priori error bounds.




35Q30,          
35Q35,          
65N12,       	
65N15,          
65N30.          




§ INTRODUCTION


Let Ø⊂ℝ^2 be an open and bounded domain with Lipschitz boundary ∂Ω. In this work, we are interested in designing and analyzing a divergence-free virtual element method (VEM) for the temperature distribution of a fluid modeled by a convection-diffusion equation coupled with a convective Brinkman–Forchheimer problem. This model can be described by the following nonlinear system of partial differential equations (PDEs):

    {[ -÷(ν(T)∇)+(·∇)
    +
    
    +
    ||^r-2
    +∇                             =                       f in Ω,;                           ÷()                             =                    0   in  Ω,;              -÷(κ(T)∇ T)+·∇ T                             =                   g 
      in  Ω, ].

together with the Dirichlet boundary conditions = 0 and T=0 on ∂Ω. The unknowns of the system are the velocity field , the
pressure , and the temperature T of the fluid. The data are the external force f, the external heat source g, the viscosity coefficient ν, and the thermal diffusion coefficient κ. We note that ν(·) and κ(·) are coefficients that can depend nonlinearly on the temperature. Finally, the parameter r is chosen such that r ∈ [3, 4]. For a discussion of this parameter r, the Forchheimer term in (<ref>), and its physical implications, we refer to the reader to <cit.> and the references therein.

Beginning with the pioneering works <cit.>, great efforts have been made to develop and analyze VEMs. These methods are a relatively new family of solution techniques that allow general polytopal meshes, arbitrary polynomial degrees, and yet conforming ^̋1-approximations. The peculiarity of VEMs is that the discrete spaces consist of functions that are not known pointwise, but about which a limited amount of information is available. This limited information is sufficient to construct stiffness matrices and right-hand sides. On the other hand, VEMs allow great flexibility with regard to the shape of the elements; for example, convex and non-convex elements are permitted. This is an advantage over the classical finite element methods (FEMs), which make it possible to deal with domains that are difficult to discretize with triangles or quadrilaterals. For these reasons, the VEM is generally considered as a generalization of the FEM. The analysis of VEMs has been successfully developed for a large number of problems. We refer the interested reader to <cit.> for a recent review and discussion. As for the development of VEMs for various linear and nonlinear flow problems, we refer the reader to the non-exhaustive list <cit.>. An important feature of some methods for the treatment of incompressible fluids is that the divergence-free condition is preserved at the discrete level.

The aim of the present work is to develop and analyze a VEM to approximate the velocity, pressure, and temperature variables that solve the system (<ref>). Our analysis is inspired by the recent work <cit.>, in which the authors propose a divergence-free VEM to solve the coupling between the Navier–Stokes equations and a suitable heat equation. Here, the viscosity coefficient depends on the temperature variable. In our work, we complement and extend the results in <cit.> in the following two directions: First, the model we consider allows the dependence of the thermal diffusion coefficient κ(·) on the temperature variable. Second, we consider the so-called convective Brinkman–Forchheimer equations. In contrast to <cit.>, our model also considers the term + ||^r-2, with r ∈ [3,4]. The consideration of this term was suggested by Forchheimer <cit.>, who realized that Darcy’s law is not adequate for moderate Reynolds numbers. Indeed, Forchheimer found that the relationship between the Darcy velocity and the pressure gradient was nonlinear and that this nonlinearity appeared to be quadratic (r=3) for a variety of experimental data <cit.>. This leads to a modification of the Darcy equations, usually referred to as the Darcy–Forchheimer equations <cit.>. In <cit.>, Forchheimer also noted that some data sets could not be described by the quadratic correction so he also postulated that the correction of Darcy’s law could allow a polynomial expression for u, e.g., u + |u|u + |u|^2u and u + |u|^r-2u; see <cit.>, <cit.>, and <cit.>. In practice, r takes the value 3 and 4 in several applications  <cit.> and also fractional values such as r = 7/2 <cit.>. This is the reason why we consider r ∈ [3,4]. In view of these considerations, (<ref>) can be seen as an extension of the model in <cit.>.



 §.§ Contributions


To the best of our knowledge, this is the first paper that analyzes a VEM for the nonlinear coupled problem (<ref>). Since several sources of nonlinearity are involved, analyzing a solution technique is far from trivial. In the following, we list what we consider to be the main contributions of our work:

∙ A VEM: Inspired by the scheme proposed in <cit.> and the discrete spaces proposed and analyzed in <cit.>, we propose the VEM (<ref>).

∙ Existence and uniqueness of discrete solutions: We derive an existence result for the discrete problem (<ref>) without restriction on the
problem data by using a fixed point strategy; see Theorem <ref>. Moreover, we obtain a global uniqueness result when the problem data is suitably restricted; see Theorem <ref>.

∙ Optimal error estimates: Assuming that the continuous and discrete problems admit unique solutions and under suitable regularity assumptions for the continuous solution, we derive optimal error estimates for the proposed VEM; see Theorem <ref>. The analysis borrows ideas and components from <cit.>, <cit.>, and <cit.>.





 §.§ Outline

The paper is structured as follows: We begin with section <ref>, where we introduce notations and basic assumptions that we will use in our work. In sections <ref>  and <ref>, we summarize some results related to the convective Brinkmann–Forchheimer problem and a suitable heat equation, respectively.


Section <ref> is devoted to review existence, uniqueness, and stability results for the coupled problem (<ref>). The core of our work begins in section <ref>, where we first introduce the standard techniques for analyzing VEMs and propose our discrete scheme (<ref>).

In section <ref>, we prove the existence of discrete solutions without restriction on the problem data, as well as a global uniqueness result when the problem data is appropriately restricted. Section <ref> is devoted to the development of a rigorous analysis of error estimates. Finally, in section  <ref>, we report a series of numerical tests in which we evaluate the performance of the proposed method for different configurations and polygonal meshes.





§ NOTATION AND PRELIMINARY REMARKS


Let us establish the notation and the framework within which we will work.

 

 §.§ Notation

 
 In this paper, Ω is an open and bounded polygonal domain of ℝ^2 with Lipschitz boundary ∂Ω. If 𝒳 and 𝒴 are normed vector spaces, we write 𝒳↪𝒴 to denote that 𝒳 is continuously embedded in 𝒴. We denote by 𝒳' and ·_𝒳 the dual and the norm of 𝒳, respectively.
For p ∈ (1,∞), we denote by q ∈ (1,∞) its Hölder conjugate, which is such that 1/p + 1/q = 1. The relation ≲ means that ≤ C, with a positive constant C that is independent of , , and the discretization parameters. The value of C can change at each occurrence. 

We use the standard notation for Lebesgue and Sobolev spaces. The spaces of vector-valued functions and the vector-valued functions themselves are denoted by bold letters. In particular, we use the following notation: :=_̋0^1(Ø), :=[_̋0^1(Ø)]^2, and :=Ł_0^2(Ø). As usual, the dual of _̋0^1(Ø) is denoted by ^̋-1(Ø). We also introduce the space
:= {∈ :  ÷ = 0}.
We conclude this section with the classical and well–known Poincaré inequality: For v ∈, there exists = (Ω) such that

    v _1,Ø≤| v |_1,Ø    ∀ v ∈.

For the sake of simplicity, we assume that the vector-valued couterpart of (<ref>), i.e., 𝐯_1,Ø≤| 𝐯 |_1,Ø for all 𝐯∈𝐕, holds with the same constant .













 §.§ Data assumptions



We make the following assumptions on the viscocity ν(·) and the diffusion coefficient κ(·).



  A0) κ(·) is extrictly positive, bounded, and Lipschitz continuous, i.e., there exist constants κ_*, κ^*, κ_lip > 0 such that 

    0 < κ_* ≤κ(r) ≤κ^*,
         
    |κ(r_1) - κ(r_2)| ≤κ_lip|r_1 - r_2| 
        ∀ r,r_1,r_2 ∈ℝ.


  A1) ν(·) is extrictly positive, bounded, and Lipschitz continuous, i.e., there exist constants ν_*, ν^*, ν_lip > 0 such that 

    0<ν_* ≤ν(r) ≤ν^*, 
        
    |ν(r_1) - ν(r_2)| ≤ν_lip|r_1 - r_2|
        ∀ r,r_1,r_2 ∈ℝ.

 




§ A CONVECTIVE BRINKMAN–FORCHHEIMER PROBLEM


We present existence and uniqueness results for the following weak formulation of a convective Brinkman–Forchheimer problem: Given f∈[^̋-1(Ø)]^2, find (,)∈× such that 

    a_L(,)+c_N(;,)+c_F(;,)+ d(,) + b(,)
    = 
    ⟨f,⟩,
       
    b(,) = 0,

for all (,)∈×. Here, ⟨· ,·⟩ denotes the duality pairing between [^̋-1(Ø)]^2 and 𝐕. The forms that occur in (<ref>) are defined as follows: a_L:×→ℝ, c_N, c_F :××→ℝ, d:×→ℝ, and b:×→ℝ are such that

    a_L(,)     := ∫_Øν∇ : ∇,
        
    c_N(;,)  :=∫_Ø (·∇)·,
    
    
    c_F(;,)     :=∫_Ø ||^r-2·,
        
    d(,) :=∫_Ø·,

and b(,):= -∫_Ø÷().

We list some of the most important properties that these forms satisfy:



  * a_L(·,·) is a coercive and continuous bilinear form: For every ,∈, we have

    a_L(,) ≥ν_*||_1,Ø^2, 
         
    a_L(,) ≤ν^*||_1,Ø||_1,Ø.



  * c_N(·;·,·) is skew-symmetric: For every ∈ and ,∈, we have

    c_N(;,) + c_N(;,)= 0,
        
    c_N(;,) = 0.

In addition, c_N(·;·,·) is continuous: For every , ,∈, we have

    c_N(;,) ≤ C_N||_1,Ø||_1,Ø||_1,Ø.



  * c_F(·;·,·) is continuous: For every , ,∈, we have

    c_F(;,) ≤ C_F||_1,Ø^r-2||_1,Ø||_1,Ø.


  * b(·,·) is a continuous bilinear form: For every ∈ and ∈, we have

    b(,) ≤ ||_1,Ø_0,Ø.

Moreover, b(·,·) satisfies an inf-sup condition: There exists a constant β>0 such that

    ∈sup b(,)||_1,Ø≥β_0,Ø  ∀∈.

 



 §.§ Existence, stability, and uniqueness results


We present existence, uniqueness, and stability results for solutions of (<ref>). In particular, we present the existence of solutions without restriction on the data and a global uniqueness result when the data is suitably restricted <cit.>.


Let Ø⊂ℝ^2 be an open and bounded domain with Lipschitz boundary ∂Ø, and let f∈[^̋-1(Ø)]^2. Then, there exists at least one solution (,)∈× for problem (<ref>) which satisfies 

    ||_1,Ø≤ν_*^-1f_-1, 
        _0,Ø≤β^-1Λ(f)f_-1,

 where Λ(f) := 1 + ν^*ν_*^-1 + C_Nν_*^-2f_-1+^2ν_*^-1 + C_Fν_*^1-rf_-1^r-2. If, in addition, C_Nf_-1 < ν_*^2, then there is a unique pair (,)∈× that solves (<ref>).







§ A NONLINEAR HEAT EQUATION


We consider the following nonlinear heat equation in weak form: Given g∈^̋-1(Ø) and ∈, find T∈ such that

    𝔞(T;T,S)+𝔠(;T,S) = ⟨ g,S⟩    ∀ S∈.

Here, ⟨·,·⟩ denotes the duality pairing between ^̋-1(Ø) and ^̋1(Ø). The forms 𝔞 and 𝔟 are defined as follows:

    𝔞:××→ℝ,     𝔞(X;R,S)   := ∫_Øκ(X)∇ R ·∇ S, 
    
    𝔠:××→ℝ,     𝔠(;R,S)   := ∫_Ø (·∇ R)S.


We list some of the most important properties that these forms satisfy:



  * Given X ∈, the form 𝔞(X;·,·) is coercive and continuous: 


    𝔞(X;S,S) ≥κ_*|S|_1,Ø^2, 
        𝔞(X;R,S) ≤κ^*|R|_1,Ø|S|_1,Ø    ∀ R,S ∈.


  * 𝔠(·;·,·) is skew-symmetric: For every ∈ and R,S ∈, we have

    𝔠(;R,S) + 𝔠(;S,R) = 0, 
        𝔠(;S,S) = 0.

In addition, 𝔠(·;·,·) is continuous: For every ∈ and R,S ∈, we have

    𝔠(;R,S) 
    ≤ℭ||_1,Ø|R|_1,Ø|S|_1,Ø.


 



 §.§ Existence, stability, and uniqueness results
 
We present existence, stability, and uniqueness results for problem (<ref>) <cit.>. 


Let Ω⊂ℝ^2 be an open and bounded domain with Lipschitz boundary ∂Ω, and let g∈^̋-1(Ø). Then, there exists at least one solution T∈_̋0^1(Ø) of problem (<ref>) which satisfies 

    |T|_1,Ø≤κ_*^-1g_-1.

If, in addition, (<ref>) has a solution T_1 ∈_𝔮^1(Ø) ∩_̋0^1(Ω) such that

    κ_*^-1κ_lipC_𝔭↪ 2|T|__𝔮^1(Ø) < 1,

then problem (<ref>) has no other solution T_2 ∈_̋0^1(Ω). Here, C_𝔭↪ 2>0 is the best constant in _̋0^1(Ø)↪Ł^𝔭(Ø), where 𝔭 < ∞, and 𝔮 is chosen such that 1/𝔭+1/𝔮=1/2.







§ THE COUPLED PROBLEM

We now introduce a weak formulation for the system (<ref>) and review existence and uniqueness results. The weak formulation is as follows: Given f∈[^̋-1(Ø)]^2 and g∈^̋-1(Ø), find (,,T)∈×× such that


    {[ a(T;,) + c_N(;,)+c_F(;,)+ d(,) + b(,)                                     =                                 ⟨f,⟩,;                                  b(,)                                     =                                    0,;                      𝔞(T;T,S)+𝔠(;T,S)                                     =                               ⟨ g,S⟩, ].

for all (,,S) ∈××. The forms c_N(·; ·,·), c_F(·; ·,·), d(·,·), and b(·,·) are defined in Section <ref> while 𝔞(·; ·,·) and 𝔠(·; ·,·) are defined in Section <ref>. The form a is defined as follows:

    a:××→ℝ,
        
    a(X;,):=∫_Øν(X)∇ : ∇.


Since ν satisfies the properties in A1), it is immediate that, given X ∈,

    a(X;,) ≥ν_*||_1,Ø^2, 
         
    a(X;,) ≤ν^*||_1,Ø||_1,Ø    ∀,∈.


Let ∈. Define the forms c_N^skew( ; ·, ·): ×→ℝ and 𝔠^skew( ; ·, ·): ×→ℝ by
c_N^skew(;,) := 1/2 [c_N(;,)-c_N(;,)]
and
𝔠^skew(;R,S) := 1/2[𝔠(;R,S)-𝔠(;S,R)], 
respectively. As discussed in Sections <ref> and <ref>, c_N and 𝔠 are skew-symmetric. As a result, c_N(;,) = c_N^skew(;,) for every ,∈ and 𝔠(;R,S) = 𝔠^skew(;R,S) for every R,S ∈. This observation will be important for the development of a virtual element numerical scheme.





 §.§ Existence, stability, and uniqueness results

We review the existence of solutions for the system (<ref>) without restriction on the data and a global uniqueness result when the data is suitably restricted and the solution is slightly smoother <cit.>.


Let Ø⊂ℝ^2 be an open and bounded domain with Lipschitz boundary ∂Ø, let f∈[^̋-1(Ø)]^2, and let g ∈^̋-1(Ω). Let ν and κ be as in Section <ref>. Then, the nonlinear system (<ref>) has at least one solution (,,T) ∈××. Moreover, we have

    ||_1,Ø≤ν_*^-1f_-1, 
      _0,Ø≤β^-1Λ(f)f_-1,
      
    |T|_1,Ø≤κ_*^-1g_-1,

where Λ(f) is defined in the statement of Proposition <ref>.  Furthermore, if the system (<ref>) has a solution (_1,_1,T_1) ∈ (∩_2+ε^1(Ø))××(∩_2+ε^1(Ø)) for some ε > 0 such that

    C_Nν_*^2f_-1 + ν_lip𝒞_ε𝒞^2_4 → 2g_-1|_1|__2+ε^1(Ø)ν_*κ_*(κ_*-κ_lipC_ϵ|T_1|__2+ε^1(Ø))   < 1,
    
    κ_lip𝒞_ϵ|T_1|__2+ε^1(Ø)   < κ_*,

then the system (<ref>) has no other solution (_2,_2,T_2)∈××. Here, 𝒞_4↪ 2 and 𝒞_ε denote the best constant in ↪Ł^4(Ø) and ↪Ł^2(2+ε)/ε(Ø) respectively.

























































































































































§ A VIRTUAL  ELEMENT APPROXIMATION


In this section, we introduce a virtual element approximation for the nonlinear system (<ref>) and derive error estimates. For this purpose, we first introduce some notions and basic ingredients <cit.>. From now on, we assume that Ω is a polygonal and bounded domain with Lipschitz boundary.

Let {_h}_h>0 be a sequence of decompositions of Ω into general polygonal elements E. Here, h:=max{ h_E: E∈_h } and h_E denotes the the diameter of the element E. We denote by |E| the area of the element E ∈_h.



 §.§ Mesh regularity

We make the following assumptions on the sequence {_h}_h>0 <cit.>: There exists ϱ > 0 such that for all h and every E in _h,


  A2) E is star-shaped with respect to a ball B_ with radius ≥ϱ h_E <cit.> and

  A3) the distance between any two vertices of E is ℓ≥ϱ h_E <cit.>.




 §.§ Basic ingredients

Let k ∈ℕ, t ∈ℝ^+, and let p ∈ [1,+∞]. Let A be an open and bounded domain in ℝ^2. We introduce some basic spaces that will be useful later. First, we introduce ℙ_k(A) — the set of polynomials on A of degree less than or equal to k —  with the convention that ℙ_-1(A) = { 0 }. Secondly, ℙ_k(𝒯_h) = {q ∈Ł^2(Ø) :  q|_E∈ℙ_k(E)  ∀ E ∈𝒯_h}. Finally, we introduce the space ^t_p(𝒯_h) = { v ∈Ł^2(Ø) :   v|_E∈^t_p(E)  ∀ E ∈𝒯_h}. We equip the space ^t_p(𝒯_h) with the broken norm ·__p^t(_h) and the broken seminorm | · |__p^t(_h), which are defined by

    [ v__p^t(_h) := [ ∑_E∈_hv__p^t(E)^p]^1/p, 
       
    |v|__p^t(_h) := [ ∑_E∈_h |v|__p^t(E)^p]^1/p, 
      
    p ∈ [1, +∞), ]

respectively. If p = + ∞, then v__∞^t(_h) :=max{v__∞^t(E): E∈_h }. We denote the vector-valued counterparts of ℙ_k(A), ℙ_k(𝒯_h), ^t_p(𝒯_h) by [ℙ_k(A)]^2, [ℙ_k(_h)]^2, and ^t_p(𝒯_h), respectively.

For an element E∈_h, we denote by x_E the centroid of E. For a multi-index α = (α_1,α_2) ∈ℕ^2, we define |α| = α_1 + α_2. Let n ∈ℕ. A natural basis associated to the space ℙ_n(E) is the set of normalized monomials

    𝕄_n(E) := {𝐦_α: α∈ℕ^2,   |α| ≤ n },
        𝐦_α(𝐱) := (𝐱-x_Eh_E)^α.

Let m ∈ℕ be such that m ≤ n. We also introduce

    ℙ_n\ m(E) := span{𝐦_α: 
    α∈ℕ^2,   m + 1 ≤ |α| ≤ n }.




 §.§ Projections

In this section we introduce some appropriate projections that will be useful for our analysis. Let E ∈_h, and let n ∈ℕ∪{0}. We introduce the Ł^2(E)-orthogonal projection as follows:

    Π^0,E_n: Ł^2(E) →ℙ_n(E):
       
     (q_n, v  - Π^0,E_nv )_Ł^2(E) = 0   ∀ q_n∈ℙ_n(E).

We denote by Π^0,E_n the vector-valued couterpart of Π^0,E_n. We also introduce the ^̋1(E)-seminorm projection, which is defined as follows:

    Π^∇,E_n: ^̋1(E) →ℙ_n(E):
       
     (∇ q_n, ∇(v  - Π^∇,E_nv) )_Ł^2(E) = 0   ∀ q_n∈ℙ_n(E)

and ∫_∂ E (v  - Π^∇,E_nv) = 0. We denote by Π^∇,E_n the vector-valued couterpart of Π^0,E_n.



 §.§ Virtual element spaces for the velocity variable


Following <cit.>, we introduce a finite-dimensional local virtual space for E ∈_h:

    (E):= {_h ∈ [^̋1(E)]^2∩ [𝒞^0(∂ E)]^2 -Δ_h - ∇∈𝐱^⊥ℙ_k-1(E), for some ∈Ł_0^2(E), 
    ÷_h ∈ℙ_k-1(E),  _h|_e∈ [ℙ_k(e)]^2 ∀ e ⊂∂ E},
        𝐱^⊥ := (x_2,-x_1).

With this space at hand, we introduce (E) as a restriction of the space (E):

    (E) := {_h ∈(E):  (_h - ^∇,E_k_h,𝐱^⊥q_k-1)_𝐋^2(E) = 0  ∀q_k-1∈ℙ_k-1\ k-3(E)}.

A similar definition of _h(E) can be found in <cit.>; compare with <cit.>.

We recall some properties of the virtual space (E) <cit.>, which are based on the description in <cit.> and <cit.>.


 
  (P1) Polynomial inclusion: [ℙ_k(E)]^2 ⊆(E);
 
  (P2) Degrees of freedom (DoFs): the following linear operators 𝐃_𝐕 constitute a set of DoFs for the virtual element space (E):
 
 
  𝐃_𝐕1 the values of 𝐯_h at the vertices of E,
 
  𝐃_𝐕2 the values of _h at k-1 distinct points of every edge e ∈∂ E,
 
  𝐃_𝐕3 the moments of _h
 
    1|E|∫_E_h·𝐦^⊥𝐦_α  ∀𝐦_α∈𝕄_k-3(E), 
        𝐦^⊥ = 𝐱^⊥-𝐱_E^⊥h_E.


  𝐃_𝐕4 the moments of ÷_h

    h_E|E|∫_E÷_h 𝐦_α  ∀𝐦_α∈𝕄_k-1(E), 
      
    |α| = α_1+α_2 > 0.

 


As a final ingredient, we introduce _k-1^0,E: ∇(E) → [ℙ_k-1(E)]^2 × 2 such that

    (q_k-1,∇_h - ^∇,E_k-1∇_h)_𝐋^2(E) = 0   ∀q_k-1∈ [ℙ_k-1(E)]^2 × 2.


The following remark is in order.

The DoFs 𝐃_𝐕 allow an exact computability of

    ^0,E_k: (E) → [ℙ_k(E)]^2,
        _k-1^0,E: ∇(E) → [ℙ_k-1(E)]^2 × 2

in the following sense: Given 𝐯_h ∈(E), ^0,E_k 𝐯_h and _k-1^0,E∇𝐯_h can be computed using only, as unique information, the DoFs values 𝐃_𝐕 of 𝐯_h <cit.>.



Finally, we define the global velocity space  as follows <cit.>:

    := {_h ∈ :  _h|_E∈(E)  ∀ E ∈_h}.




 §.§ Finite element space for the pressure variable


We introduce the following finite element space for approximating a pressure variable:

    _h := {_h ∈ _h|_E∈ℙ_k-1(E)  ∀ E ∈_h}.

It is important to note that the pair (,_h) satisfies the following discrete inf-sup condition: there exists a constant β̃>0 independent of h
such that

    _h ∈sup b(_h,_h)|_h|_1,Ø≥β̃_h_0,Ø  ∀_h ∈_h;

<cit.>. Let us now introduce the discrete kernel

    := {_h ∈:   b(_h,_h) = 0   ∀_h ∈_h}.


The following observation is important:

Let _h ∈. Given that ÷_h|_E ∈ℙ_k-1(E), we deduce that 𝐙_h ⊆𝐙: the functions in the discrete kernel are exactly divergence-free.






 §.§ Virtual element spaces for the temperature variable


To introduce a virtual element space for approximating a temperature variable, we first introduce, for each E ∈_h, the finite-dimensional local virtual space <cit.>

    (E):= {S_h ∈^̋1(E) ∩𝒞^0(∂ E) Δ S_h ∈ℙ_k(E),  S_h|_e∈ℙ_k(e)  ∀ e ∈∂ E}.

With this space at hand, on each E ∈_h we introduce <cit.>

    (E) := {S_h ∈(E):   (S_h - Π^∇,E_k S_h, q_k)_Ł^2(E) = 0   ∀ q_k ∈ℙ_k∖ k-2(E) }.



We recall some properties of the space (E) based on the presentation in <cit.> and <cit.>:



  (P3) Polynomial inclusion: ℙ_k(E) ⊆(E);

  (P4) Degrees of freedom: the following linear operators D_V constitute a set of DoFs for the virtual element space (E):


  D_V1 the values of S_h at the vertices of E,

  D_V2 the values of S_h at k-1 distinct points of every edge e ∈∂ E,

  D_V3 the moments of S_h

    1|E|∫_E S_h m_α  ∀m_α∈𝕄_k-2(E).




The DoFs D_V allow an exact computability of 

    Π^0,E_k: _h(E) →ℙ_k(E),
         ^0,E_k-1: ∇_h(E) → [ℙ_k-1(E)]^2,

<cit.>; compare with Remark <ref>.


 
Finally, we introduce a global virtual space to approximate a temperature variable

    := {S_h ∈:   S_h|_E∈(E)  ∀ E ∈_h}.




 §.§ Virtual element forms


With the discrete spaces , _h, and   in hand, and following <cit.> we now introduce discrete versions of the continuous forms involved in the weak problem (<ref>).




  * The discrete counterpart of a(·;·,·) (cf. (<ref>)) is defined by 

    a_h(·;·,·):××→ℝ, 
        
    a_h(X_h;_h,_h) := ∑_E∈_h a_h^E(X_h;_h,_h),

where a_h^E(·;·,·):(E)×(E)×(E) →ℝ is given by

    a_h^E(X_h;_h,_h) := 
    ∫_E ν(Π_k^0,EX_h)^0,E_k-1∇_h:^0,E_k-1∇_h 
    
    
    + 
    ν(Π^0,E_0X_h)S_V^E((𝐈-^0,E_k)_h,(𝐈-^0,E_k)_h).

Here, S_V^E(·,·):(E)×(E) →ℝ is a computable symmetric form that satisfies

    |_h|_1,E^2 ≲ S_V^E(_h,_h) ≲ |_h|_1,E^2   ∀_h ∈(E) ∩(^0,E_k).


  * The discrete counterpart of 𝔞(·;·,·) (cf. (<ref>)) is defined by

    𝔞_h(·;·,·):××→ℝ, 
        𝔞_h(X_h;R_h,S_h) := ∑_E∈_h𝔞_h^E(X_h;R_h,S_h),

where 𝔞_h^E(·;·,·):(E)×(E)×(E) →ℝ is given by

    𝔞_h^E(X_h;R_h,S_h) := ∫_E κ(Π^0,E_k X_h)^0,E_k-1∇ R_h·^0,E_k-1∇ S_h 
    
    + κ(Π^0,E_0 X_h)S_T^E((-Π^0,E_k)R_h,(-Π^0,E_k)S_h).

Here, S_T^E(·,·):(E)×(E) →ℝ is a computable symmetric form that satisfies

    |X_h|_1,E^2 ≲ S_T^E(X_h,X_h) ≲ |X_h|_1,E^2   ∀ X_h ∈(E) ∩(Π^0,E_k).



  * The discrete counterpart of c_N(·;·,·) (cf. (<ref>)) is defined by

    c_N,h(·;·,·):××→ℝ, 
         
    c_N,h(_h;_h,_h) := ∑_E∈_h c_N,h^E(_h;_h,_h),

where c_N,h^E(·;·,·):(E)×(E)×(E) →ℝ is given by

    c_N,h^E(_h;_h,_h) := 
    ∫_E [^0,E_k-1(∇_h) ^0,E_k _h]·^0,E_k _h.

Following <cit.> and <cit.>, we define

    c_N,h^skew,E(_h;_h,_h):=12(c_N,h^E(_h;_h,_h)-c_N,h^E(_h;_h,_h)),

and c_N,h^skew(_h;_h,_h) := ∑_E∈_h c_N,h^skew,E(_h;_h,_h).


  * The discrete counterpart of c_F(·;·,·) (cf. (<ref>)) is defined by

    c_F,h(·;·,·):××→ℝ, 
         
    c_F,h(_h;_h,_h) := ∑_E∈_h c_F,h^E(_h;_h,_h),

where c_F,h^E(·;·,·):(E)×(E)×(E) →ℝ is given by 

    c_F,h^E(_h;_h,_h) := c_F^E(^0,E_k_h;^0,E_k_h,^0,E_k_h).



  * The discrete counterpart of 𝔠(·;·,·) (cf. (<ref>)) is defined by

    𝔠_h(·;·,·):××→ℝ, 
        𝔠_h(_h;R_h,S_h) := ∑_E∈_h𝔠_h^E(_h;R_h,S_h),

where 𝔠_h^E(·;·,·):(E)×(E)×(E) →ℝ is given by

    𝔠_h^E(_h;R_h,S_h) := ∫_E (^0,E_k_h ·^0,E_k-1∇ R_h)Π^0,E_kS_h.

Following <cit.>, we define 

    𝔠_h^skew,E(_h;R_h,S_h):=12(𝔠_h^E(_h;R_h,S_h)-𝔠_h^E(_h;S_h,R_h)),

and 𝔠_h^skew(_h;R_h,S_h) := ∑_E∈_h𝔠_h^skew,E(_h;R_h,S_h).

  * The discrete counterpart of d(·,·) (cf. (<ref>)) is defined by

    d_h(·,·) : ×→ℝ, 
         d_h(_h,_h) := ∑_E∈_h d_h^E(_h,_h),

where d_h^E(·,·):(E)×(E) →ℝ is given by

    d_h^E(_h,_h) := ∫_E ^0,E_k_h ·^0,E_k_h.



We note that c_N,h^skew and 𝔠_h^skew are skew-symmetric by construction, i.e., for _h ∈_h, 

    c_N,h^skew,E(_h;_h,_h) = 0,
         𝔠_h^skew(_h;R_h,R_h) = 0,

for every _h ∈_h and every R_h ∈_h. As usual, this property simplifies the analysis of the corresponding discrete problem <cit.>.


As mentioned in <cit.>, we do not introduce any approximation of the bilinear form b. We note that b(_h,_h) for _h ∈_h and _h∈_h is computable from the DoFs 𝐃_𝐕1, 𝐃_𝐕2, 𝐃_𝐕4 because _h is a polynomial on each element E ∈_h.




 §.§ Virtual element forcing


From now on we will assume that f∈𝐋^2(Ω) and g ∈Ł^2(Ω). Within this framework, we define the discrete sources f_h ∈ [ℙ_k(_h)]^2 and g_h ∈ℙ_k(_h) to be such that, for each E ∈_h,

    f_h|_E := ^0,E_kf, 
         
    g_h|_E := Π^0,E_k g.




 §.§ The virtual element method

With all the ingredients and definitions introduced in the previous sections, we finally design a virtual element discretization for the nonlinear system (<ref>): Find  (_h,_h,T_h)∈×_h× such that


    {[ a_h(T_h;_h,_h)
    +
    c_N,h^skew(_h;_h,_h)
    +
    c_F,h(_h;_h,_h);                                 +
    d_h(_h,_h) + b(_h,_h)                                                       =                                           (f_h,_h)_0,Ø,;                                                b(_h,_h)                                                       =                                                      0,;                 𝔞_h(T_h;T_h,S_h)
    +
    𝔠_h^skew(_h;T_h,S_h)                                                       =                                          (g_h,S_h)_0,Ø, ].

for all (_h,_h,S_h) ∈×_h×.


Using the space  defined in (<ref>), the problem (<ref>) can be reformulated as follows: Find  (_h,T_h)∈× such that

    {[ a_h(T_h;_h,_h)+c_N,h^skew(_h;_h,_h)+c_F,h(_h;_h,_h);                                        + d_h(_h,_h)                                                   =                                       (f_h,_h)_0,Ø,;               𝔞_h(T_h;T_h,S_h)+𝔠_h^skew(_h;T_h,S_h)                                                   =                                      (g_h,S_h)_0,Ø, ].

for all (_h,S_h) ∈×.


To provide an analysis for the discrete problem (<ref>), we list some properties that the discrete forms satisfy:



  * For any X_h ∈_h, a_h(X_h;·,·) is a coercive and continuous bilinear form: There exist α_*>0 and α^*>0 such that, for every _h,_h ∈ we have

    a_h(X_h;_h,_h) ≤α^*ν^*|_h|_1,Ø|_h|_1,Ø,
        
    a_h(X_h;_h,_h) ≥α_*ν_*|_h|_1,Ø^2.

These inequalities result from standard properties of Π^0,E_n and _k-1^0,E in conjunction with basic inequalities and the properties A1 and (<ref>) that ν and S_V satisfy.


  * For any X_h ∈_h, 𝔞_h(X_h;·,·) is a coercive and continuous bilinear form: There exist β_*>0 and β^*>0 such that, for every R_h, S_h ∈, we have

    𝔞_h(X_h;R_h,S_h) ≤β^*κ^*|R_h|_1,Ø|S_h|_1,Ø,
        𝔞_h(X_h;S_h,S_h) ≥β_*κ_*|S_h|_1,Ø^2.

These inequalities follow from very similar arguments to those that lead to (<ref>).


  * c_N,h^skew(·;·,·) is continuous: For every _h,_h,_h ∈, we have

    c_N,h^skew(_h;_h,_h) ≤C_N|_h|_1,Ø|_h|_1,Ø|_h|_1,Ø.

A proof of this estimate is essentially contained in <cit.>.


  * c_F,h(·;·,·) is continuous: For every _h,_h,_h ∈, we have

    c_F,h(_h;_h,_h) ≤C_F|_h|_1,Ø^r-2|_h|_1,Ø|_h|_1,Ø.

A proof of this estimate can be found in <cit.>.


  * 𝔠_h^skew(·;·,·) is continuous: For every _h ∈_h and R_h,S_h ∈_h, we have

    𝔠_h^skew(_h;R_h,S_h) ≤ℭ|_h|_1,Ø |R_h |_1,Ø|S_h|_1,Ø,

This bound can be obtained with arguments similar to those in the proof of <cit.>.


  * d_h(·,·) is continuous: For every _h, _h ∈_h, we have

    d_h(_h,_h) ≤_h_0,Ø_h_0,Ø.

This bound follows directly from the Cauchy-Schwarz inequality in 𝐋^2(E), the 𝐋^2-stability of _k^0,E, and the Cauchy-Schwarz inequality in ℝ^#_h.


  * b(·,·) is continuous: For every _h ∈_h and _h ∈_h, we have

    b(_h,_h) ≤ |_h|_1,Ø_h_0,Ø.
 

 



§ EXISTENCE AND UNIQUENESS OF SOLUTIONS FOR OUR VEM


In this section, we derive an existence result for the discrete problem (<ref>) without restriction on the problem data by using a fixed point strategy. Moreover, we obtain a global uniqueness result when the problem data is suitably restricted.


Under the data assumptions A0) and A1), the nonlinear system (<ref>) admits at least one solution (_h,_h,T_h) ∈×_h×. Moreover, the following estimates hold uniformly in h:

    |_h|_1,Ø   ≤α_*^-1ν_*^-1f_h_0,Ø,
      _h_0,Ø≤β̃^-1Γf_h_0,Ø, 
    
    |T_h|_1,Ø   ≤β_*^-1κ_*^-1g_h_0,Ø,

Here,
Γ :=
1
+
(α_*ν_*)^-1(
α^*ν^*
+
C_Nα_*^-1ν_*^-1f_h_0,Ø
+
^r-2C_Fα_*^2-rν_*^2-rf_h_0,Ø^r-2
+
^2).


 We proceed in two steps.

Step 1. Existence of solutions: Let us first analyze the existence of solutions to the reduced problem (<ref>). To this purpose, we

define A_(_h,T_h): ×→ℝ, for (_h,T_h) ∈×, as follows:

    A_(_h,T_h)(_h,S_h) := a_h(T_h;_h,_h)+c_N,h^skew(_h;_h,_h) + c_F,h(_h;_h,_h) + d_h(_h,_h)
    
    
    - (f_h,_h)_0,Ω + 𝔞_h(T_h;T_h,S_h) + 𝔠_h^skew(_h;T_h,S_h) - (g_h,S_h)_0,Ω.

Set (_h,S_h) = (_h,T_h) and use the skew-symmetry of c_N,h^skew(·;·,·) and 𝔠_h^skew(·;·,·) and the coercivity properties in (<ref>) and (<ref>) to obtain

    A_(_h,T_h)(_h,T_h) ≥α_*ν_*|_h|_1,Ø^2 + β_*κ_*|T_h|_1,Ø^2 + c_F,h(_h;_h,_h) + d_h(_h,_h) 
    
    - f_h_0,Ø_h_0,Ø-g_h_0,ØT_h_0,Ø,

Since d_h(_h,_h)≥ 0 and c_F,h(_h;_h,_h) ≥ 0, we thus obtain that 

    A_(_h,T_h)(_h,T_h) ≥α_*ν_*|_h|_1,Ø^2 + β_*κ_*|T_h|_1,Ø^2 - ( f_h_0,Ø|_h|_1,Ø - g_h_0,Ø|T_h|_1,Ø)
    
    ≥ 
    (min{α_*,β_*}min{ν_*,κ_*})
    |(_h,T_h)|_1,Ø^2 - (f_h,g_h)_0,Ø
    |(_h,T_h)|_1,Ø,

where |(_h,S_h)|_1,Ø := (|_h|_1,Ø^2 + |S_h|_1,Ø^2)^1/2 and (f,g)_0,Ø := (f_0,Ø^2 + g_0,Ø^2)^1/2 for every (_h,S_h) ∈_h ×_h and (f,g) ∈𝐋^2(Ω) ×Ł^2(Ω). Hence, we obtain that A_(_h,T_h)(_h,T_h) ≥ 0 for every (_h,T_h) ∈× such that

    |(_h,T_h)|_1,Ø ^2 =^2(f_h,g_h)_0,Ø^2(min{α_*,β_*}min{ν_*,κ_*})^2 =: μ^2.
 
With reference to <cit.>, there thus exists (_h,T_h) ∈_h ×_h such that A_(_h,T_h)(_h,S_h) = 0 for all (_h,S_h) ∈_h ×_h and |(_h,T_h)|_1,Ø≤μ≲(f_h,g_h)_0,Ω; i.e., the nonlinear system (<ref>) has at least one solution (_h,T_h)∈×. Finally, the existence of a solution for the system (<ref>) follows immediately from the inf-sup condition (<ref>).

Step 2. Stability bounds: Let us first derive the bound for the velocity field _h in (<ref>). To to this, we set (_h,_h)=(_h,0) in (<ref>) and use the coercivity property in (<ref>), the skew-symmetry of c_N,h^skew(·;·,·), d(_h,_h)≥ 0, and c_F,h(_h;_h,_h)≥ 0 to obtain the following bounds:

    α_*ν_*|_h|_1,Ø^2 ≤ a_h(T_h;_h,_h)+c_F,h(_h;_h,_h)+d_h(_h,_h) 
    
    = ( f_h,_h)_0,Ø≤f_h _0,Ω |_h|_1,Ø.

This immediately yields the the desire bound for the velocity field. The estimate for the temperature follows similar arguments. Finally, the estimate for the pressure can be obtained using (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>). This concludes the proof.


In the following, we present two elementary results that are useful to show the uniqueness of the solutions of system (<ref>).




Let X_h, S_h ∈, and let _h ∈ such that _h ∈_q^1(_h) for some q>2. Then, there exists C_V>0 depending on Ø, the polynomial degree k, and the shape regularity constant ϱ such that, for any _h ∈,

    | a_h(X_h;_h,_h) - a_h(S_h;_h,_h)|≤ C_Vν_lip|X_h-S_h|_1,Ø|_h|__q^1(_h)|_h|_1,Ø.




See <cit.>.
































Let X_h, S_h ∈, and let R_h ∈ such that R_h ∈_q^1(_h) for some q>2. Then, there exists C_T>0 depending on Ø, the polynomial degree k, and the shape regularity constant ϱ such that, for any P_h ∈,

    |𝔞_h(X_h;R_h,P_h) - 𝔞_h(S_h;R_h,P_h)|≤ C_Tκ_lip|X_h-S_h|_1,Ø|R_h|__q^1(_h)|P_h|_1,Ø.



The proof is similar to that of Lemma <ref>, so we skip the details.


We are now in a position to show the uniqueness of solutions for system (<ref>).

 
Let us assume that assumptions A0) and A1) hold. If

    C_Vν_lipℭg_h_0,Ø |_h|__q^1(_h)α_*ν_*β_*κ_*(β_*κ_*-C_Tκ_lip|T_h|__q^1(_h))
    +
    C_Nα_*^2ν_*^2f_h_0,Ø    < 1,
    
    
    C_Tκ_lip|T_h|__q^1(_h)   < β_*κ_*,

and the nonlinear system (<ref>) admits a solution (_h,_h,T_h) ∈×_h ×_h such that _h ∈_q^1(_h) and T_h ∈_q^1(_h) for some q>2, then this solution is unique.




We begin the proof with the assumption that there is another solution (_h,_h,T_h) ∈×_h× of system (<ref>). Define

    _h := _h-_h ∈,
        T_h := T_h-T_h ∈,
        _h := _h-_h ∈_h.



Step 1. A bound for |T_h|_1,Ø. We begin by noting that it follows directly from the definition of 𝔠_h^skew(·;·,·) that

    𝔠_h^skew(_h;T_h,T_h) = 0 𝔠_h^skew(_h;T_h,T_h) = 𝔠_h^skew(_h;T_h,T_h).

We now use the fact that (_h,T_h) and (_h,T_h) solve the system (<ref>) to set _h = 0 and S_h = T_h ∈_h in the corresponding systems and obtain

    𝔞_h(T_h;T_h,T_h) + 𝔠_h^skew(_h;T_h,T_h) = 𝔞_h(T_h;T_h,T_h) + 𝔠_h^skew(_h;T_h,T_h).

We add and subtract the term 𝔞_h(T_h;T_h,T_h) and use (<ref>) to obtain

    𝔞_h(T_h;T_h,T_h) = [𝔞_h(T_h;T_h,T_h)-𝔞_h(T_h;T_h,T_h)] - 𝔠_h^skew(_h;T_h,T_h).

We then use the coercivity property in (<ref>), the bound (<ref>), and the estimate of Lemma <ref> to obtain

    β_*κ_*|T_h|_1,Ø^2≤ C_Tκ_lip|T_h|_1,Ø^2|T_h|__q^1(_h) + ℭ|_h|_1,Ø|T_h|_1,Ø|T_h|_1,Ø.

This yields
(β_*κ_* - C_Tκ_lip|T_h|__q^1(_h))|T_h|_1,Ø≤ℭ|_h|_1,Ø|T_h|_1,Ø.
Note that β_*κ_* - C_Tκ_lip|T_h|__q^1(_h)>0. We now invoke the stability bound for T_h in (<ref>) to obtain

    |T_h|_1,Ø≤ℭg_h_0,Ωβ_*κ_*(β_*κ_* - C_Tκ_lip|T_h|__q^1(_h))|_h|_1,Ø.


Step 2. We now obtain a bound for _h. To do this, we first note that, due to the skew-symmetry of c_N,h^skew(·;·,·), we have the following

    c_N,h^skew(_h;_h,_h) = 0 
     c_N,h^skew(_h;_h,_h) = c_N,h^skew(_h;_h,_h).

We now use the fact that (_h,T_h) and (_h,T_h) solve system (<ref>) to set _h = _h ∈_h in the first equation of the corresponding systems and obtain

    a_h(T_h;_h,_h) + c_N,h^skew(_h;_h,_h) + c_F,h(_h;_h,_h) + d_h(_h,_h) 
    
    = a_h(T_h;_h,_h) + c_N,h^skew(_h;_h,_h) + c_F,h(_h;_h,_h) + d_h(_h,_h).

We add and subtract the term a_h(T_h;_h,_h) and use (<ref>) to obtain

    a_h(T_h;_h,_h) + [c_F,h(_h;_h,_h) - c_F,h(_h;_h,_h)] + d_h(_h,_h)
    
     = [a_h(T_h;_h,_h) - a_h(T_h;_h,_h)] - c_N,h^skew(_h;_h,_h).

We then use the left-hand side bound in (<ref>), <cit.>, the bound in Lemma <ref>, (<ref>), and the fact that d_h(_h,_h)≥ 0 to obtain

    α_*ν_*|_h|_1,Ø^2≤ C_Vν_lip|T_h|_1,Ø|_h|__q^1(_h)|_h|_1,Ø + C_N|_h|_1,Ø^2|_h|_1,Ø.

We now use the bound for _h in Theorem <ref> and (<ref>) to arrive at

    α_*ν_*|_h|^2_1,Ø   ≤C_Vν_lipℭg_h_0,Ø|_h|__q^1(_h)β_*κ_*(β_*κ_*-C_Tκ_lip|T_h|__q^1(_h))|_h|^2_1,Ø + C_Nα_*ν_*f_h_0,Ø|_h|^2_1,Ø.

The previous inequality allows us to conclude that 

    (1-C_Vν_lipℭg_h_0,Ø |_h|__q^1(_h)α_*ν_*β_*κ_*(β_*κ_*-C_Tκ_lip|T_h|__q^1(_h))
    -
    C_Nα_*^2ν_*^2f_h_0,Ø)|_h|^2_1,Ø≤ 0.


Step 3. In view of (<ref>), we immediately conclude that _h = 0. We now invoke (<ref>) to obtain that T_h = 0. Finally, the inf-sup condition (<ref>) and the arguments developed in Step 2 show that _h = 0. This concludes the proof.












§ ERROR ESTIMATES FOR OUR VEM


In the following, we derive error bounds for the formulation with virtual elements (<ref>). For this purpose, we will make the following regularity assumptions:



  A4) The solutions (,,T) ∈×× and the data f,g,κ,ν of system (<ref>) satisfy the following regularity properties for some 0<s≤ k:


  i) ∈𝐇^s+1(Ø), ∈^̋s(Ø), and T∈^̋s+1(Ø).

  ii) f∈𝐇^s+1(Ø) and g ∈^̋s+1(Ø).

  iii) ν(T) and κ(T) belong to ^s_∞(Ø).



Before deriving a priori error estimates, we recall some preliminary approximation properties.





Assume that A2) and A3) hold. Let ∈∩𝐇^s+1(Ø). Then, there exists _I∈ such that

    - _I_0,E + h| - _I|_1,E≲ h_E^s+1||_s+1,D(E).

for all E ∈_h. Here, D(E) denote the union of the polygons in _h intersecting E and 0<s≤ k. Moreover, if ∈, then _I ∈.


A proof of the error bound follows from the arguments in the proof of <cit.> in combination with the error bound of <cit.>.




Assume that A2) and A3) hold. Let S ∈∩^̋s+1(Ø). Then, there exists S_I∈ such that

    S - S_I_0,E + h|S - S_I|_1,E≲ h_E^s+1|S|_s+1,D(E).

for all E ∈_h. Here, D(E) denote the union of the polygons in _h intersecting E and 0<s≤ k.


A proof of this error bound can be found in <cit.>.


We will also make use of a Bramble–Hilbert Lemma <cit.>: Let 0 ≤ t ≤ s ≤ℓ+1 and 1 ≤ q,p ≤∞ such that s-2/p > t-2/q. Then, for any E∈_h,

    q_ℓ∈ℙ_ℓ(E) inf |S - q_ℓ |_W_q^t(E)≲ h_E^s-t+2/q-2/p|S|_W_p^s(E)    ∀ S ∈ W_p^s(E).













To simplify the presentation of the material, we will write the continuous and discrete stability bounds for (,T) and (_h,T_h) as in <cit.>. This is

    ||^2_1,Ø+|T|^2_1,Ø≤ C^2_estC^-2_data,
        
    |_h|^2_1,Ø+|T_h|^2_1,Ø≤C^2_estC^-2_data.

As in <cit.>,  we will also use the following notation here: We denote by 𝔄(·), 𝔅(·), ℭ(·), 𝔇(·), etc. generic constants that are independent of h, but may depend on appropriate norms of , p, T, ν, κ, f, g, Ω, or the polynomial degree k.

We are now ready to state and prove the main result of this section.


Let us assume that the assumptions A0), A1), A2), A3), and A4) hold. Let us also assume that the smallness assumptions of Proposition <ref> and Theorem <ref> hold. Let (,,T) ∈×× and (_h,_h,T_h) ∈_h×_h×_h be the unique solutions of problems (<ref>) and (<ref>), respectively. If

    C_final,T := CC_est/C_data(
    β_*κ_* - C_lip(T)/C_sol(T)
    )^-1 > 0,
    
    C_final,^-1:= ( α_*ν_* - C_N C_est/C_data - C_for -  ^2 - C_lip()/C_sol()C_final,T) > 0,

then the following a priori error estimates hold

    |-_h|_1,Ø+|T-T_h|_1,Ø   ≤𝔐(ν,κ,,T)h^s
    +
    𝔑(f,g)h^s+2
    - _h_0,Ø   ≤𝔒(ν,κ,,T,p)h^s
    +
    𝔑(f,g)h^s+2.





We follow the proof of <cit.> and proceed in several steps.

Step 1. Interpolation error estimates.
We introduce _I ∈_h and T_I ∈_h as the interpolants of ∈ and T ∈ given by Lemmas <ref> and <ref>, respectively. We also introduce _I ∈_h as follows: for each E ∈_h, _I|_E := Π^0,E_k-1. As a direct consequence of Lemmas <ref> and <ref> and the Bramble-Hilbert bound (<ref>), we obtain

    |-_I|_1,Ø≲ h^s||_s+1,Ø,    |T-T_I|_1,Ø≲ h^s|T|_s+1,Ø,   - _I_0,Ø≲ h^s| |_s,Ø.


Define e_h := _I-_h, E_h := T_I-T_h, and _h :=  _I- _h. Since ∈, we have that _I ∈_h and thus that e_h ∈_h. In the following, we bound e_h, E_h, and _h.

Step 2. An estimate for E_h. We start with the coercivity bound in (<ref>), the definition of E_h, namely E_h := T_I-T_h, and add and subtract 𝔞(T;T,E_h) to obtain

    β_*κ_*|E_h|_1,Ø^2   ≤𝔞_h(T_h;E_h,E_h) = 𝔞_h(T_h;T_I,E_h) - 𝔞_h(T_h;T_h,E_h)
    
        = 𝔞_h(T_h;T_I,E_h) - 𝔞(T;T,E_h) + 𝔞(T;T,E_h) - 𝔞_h(T_h;T_h,E_h).

We now use the third equations of the continuous and discrete systems, (<ref>) and (<ref>), respectively, to arrive at the following estimate:

    β_*κ_*|E_h|_1,Ø^2≤[𝔞_h(T_h;T_I,E_h) - 𝔞(T;T,E_h)]
    
    
    +
    [ 𝔠_h^skew(_h;T_h,E_h) - 𝔠(;T,E_h)]
    + ( g-g_h,E_h)_0,Ω =:  +  + .


Step 2.1. We estimate . To do so, we first rewrite  using the definitions of 𝔞 and 𝔞_h, given in (<ref>) and (<ref>), respectively, and a localization argument and obtain

    =
    ∑_E∈_h{∫_Eκ(Π^0,E_k T_h)
    (^0,E_k-1∇ T_I - ∇ T)
    ·^0,E_k-1∇ E_h
    -  ∫_Eκ(T)∇ T·∇ E_h
    
    
    + ∫_Eκ(Π^0,E_kT_h)∇ T ·^0,E_k-1∇ E_h
    +
    κ(Π^0,E_0 T_h) S_T^E((-Π^0,E_k)T_I,(-Π^0,E_k) E_h)}.

We construct further differences as follows:

    = ∑_E∈_h{∫_Eκ(Π^0,E_k T_h)(^0,E_k-1∇ T_I - ∇ T)·^0,E_k-1∇ E_h
    
          + ∫_E
    (
    κ(Π^0,E_k T_h) - κ(T)
    )∇ T
    ·^0,E_k-1∇ E_h
    -
    ∫_E κ(T)∇ T · ( ∇ E_h - ^0,E_k-1∇ E_h)
    
         
    +
    κ(Π^0,E_0 T_h) S_T^E((-Π^0,E_k)T_I,(-Π^0,E_k) E_h)
    }
    =: ∑_E∈_h ( _1^E + _2^E - _3^E + _4^E).

Using the definition of ^0,E_k-1, _3^E can be rewritten as _3^E = ∫_E (𝐈 - ^0,E_k-1)(κ(T)∇ T) ·∇ E_h. As in <cit.>, the terms _1^E, _3^E, and _4^E can be controlled simultaneously under the assumption A0) and (<ref>). In fact, we have

    ∑_E∈_h  (_1^E - _3^E + _4^E)
     ≤ C
     ∑_E∈_h( κ^*^0,E_k-1∇ T_I - ∇ T _0,E.
     
    .
     +
      (𝐈 - ^0,E_k-1)(κ(T)∇ T) _0,E
     +
     κ^*∇ ( - Π^0,E_k)T_I _0,E)
     ∇ E_h _0,E,

where we used ^0,E_k-1∇ E_h _0,E≤∇ E_h _0,E and ∇ (   - Π^0,E_k) E_h _0,E≲∇ E_h _0,E. We now use the triangle inequality, Lemma <ref>, and a basic error bound for the 𝐋^2(E)-projection (see, e.g. <cit.>) to obtain

    ^0,E_k-1∇ T_I - ∇ T _0,E≤∇ (T - T_I) _0,E
     +
      (𝐈 - ^0,E_k-1) ∇ T _0,E≲ h_E^s |T|_1+s,D(E).

An estimate for ∇ ( - Π^0,E_k)T_I _0,E can be derived in view of similar arguments:

    ∇ ( - Π^0,E_k)T_I _0,E≤∇ ( - Π^0,E_k) (T_I-T)  _0,E
     +
     ∇ ( - Π^0,E_k) T _0,E
    ≲∇ (T_I-T)  _0,E
     +
    ∇ ( - Π^0,E_k) T _0,E≲ h_E^s |T|_1+s,D(E).

Finally, we bound (𝐈 - ^0,E_k-1)(κ(T)∇ T) _0,E as follows:


    (𝐈 - ^0,E_k-1)(κ(T)∇ T) _0,E≲
    h_E^s | κ(T)∇ T |_s,E≲
    h_E^s κ(T) _W_∞^s(E) | T |_s+1,E.

If we substitute the estimates (<ref>), (<ref>), and (<ref>) into (<ref>) we arrive at the bound

    ∑_E∈_h  (_1^E - _3^E + _4^E)
     ≤𝔄(κ,T) h^s ∇ E_h _0,Ω,

 where 𝔄 = 𝔄(κ,T) is a constant that depends on κ and T.

 We now turn to the derivation of a bound for _2^E. For this purpose, we invoke Hölder's inequality (1/p + 1/q +1/2 =1, where q = 2 + ε >2 is given by the uniqueness assumptions of Proposition <ref>) and the stability of the Ł^2(E)-projection in Ł^2(E) and Ł^p(E) to obtain
 
    _2^E
    ≤κ_lipT-Π^0,E_kT_h_Ł^p(E)|T|__q^1(E)∇ E_h_0,E≤κ_lip ( T- Π^0,E_kT _Ł^p(E)
    
    + Π^0,E_k(T - T_h) _Ł^p(E))|T|__q^1(E)∇ E_h_0,E≲κ_lip
    ( T- Π^0,E_kT _Ł^p(E)
    
    +  T - T_h _Ł^p(E))|T|__q^1(E)∇ E_h_0,E.

A standard error bound yields T- Π^0,E_kT _Ł^p(E)≲ h_E^s |T|__s^p(E). On the other hand, we have T - T_h _Ł^p(E)≤ T - T_I _Ł^p(E) +  T_I - T_h _Ł^p(E). Substituting these estimates into (<ref>), we obtain

    _2^E ≲κ_lip
    (
    h_E^s |T|__s^p(E) +  T - T_I _Ł^p(E) +  E_h _Ł^p(E)
    )|T|__q^1(E)∇ E_h_0,E,

where we also used that E_h = T_I - T_h. If we sum over all the elements in E in _h and apply Hölder's inequality (1/p + 1/q +1/2 =1), we obtain

    _2 := ∑_E∈_h_2^E ≲κ_lip
    (
    h^s |T|__s^p(Ω) +  h^s| T |_s+1,Ω+  E_h _Ł^p(Ω)
    )|T|__q^1(Ω)∇ E_h_0,Ω

where we have used that T - T_I _Ł^p(Ω)≲∇(T - T_I) _0,Ω≲ h^s  | T |_s+1,Ω, which follows from a basic Sobolev embedding and the error bound in Lemma <ref>. As a result,

    _2 ≤C_lip(T)/C_sol(T)∇ E_h^2_0,Ω + 𝔄(κ,T)h^s∇ E_h_0,Ω.

C_sol(T) comes from Prop. <ref>: |T|__q^1(Ω) < C_sol(T)^-1 := (κ_lip𝒞_q/κ_*)^-1,

and C_lip(T)>0 is a suitable constant. In a final step, we combine (<ref>) with (<ref>) and obtain

    =
    ∑_E∈_h  (_1^E + _2^E - _3^E + _4^E) ≤C_lip(T)/C_sol(T)∇ E_h^2_0,Ω + 𝔄(κ,T)h^s∇ E_h_0,Ω.


Step 2.2. Let us now control  in (<ref>). As a first step, we note that according to Remark <ref>, 𝔠(,T;E_h) = 𝔠^skew(,T;E_h) because ∈. We use this property, the fact that 𝔠_h^skew(_h;T_h,E_h) = 𝔠_h^skew(_h;T_h -T_I,E_h) + 𝔠_h^skew(_h;T_I,E_h) = 𝔠_h^skew(_h;T_I,E_h), and add subtract 𝔠_h^skew(;T,E_h) to rewrite  as follows:

    = 𝔠_h^skew(_h;T_I - T,E_h)
     +
     𝔠_h^skew(_h - ;T,E_h)
     +
     𝔠_h^skew(;T,E_h)
     -
     𝔠^skew(;T,E_h).

Define _1:= 𝔠_h^skew(_h;T_I - T,E_h)
 +
 𝔠_h^skew(_h - ;T,E_h). To bound _1, we first use the estimate (<ref>) and obtain

    II_1
    ≤ℭ(
    |_h|_1,Ø|T-T_I|_1,Ø + |-_h|_1,Ø|T|_1,Ø)
    |E_h|_1,Ø.

We now use the interpolation error bounds of Lemmas <ref> and <ref> to arrive at

    II_1
    ≤ℭ C_est/C_data
    |e_h|_1,Ø|E_h|_1,Ø
    +
    𝔅(u,T)h^s|E_h|_1,Ø,

where we also used (<ref>). Define _2^E:= 𝔠_h^skew,E(;T,E_h)-𝔠^skew,E(;T,E_h), for E ∈_h, and _2:= 𝔠_h^skew(;T,E_h)
 - 𝔠^skew(;T,E_h). We note that

    _2^E
    =
    1/2[
    ∫_E ( ^0,E_k ·^0,E_k-1∇ T ) Π^0,E_k E_h
    -
    ∫_E  ( ·∇ T ) E_h
    ]
    
    
    -
    1/2[
    ∫_E ( ^0,E_k ·^0,E_k-1∇ E_h ) Π^0,E_k T
    -
    ∫_E  ( ·∇ E_h ) T
    ]
    =:
    1/2_2,a^E - 1/2_2,b^E.

The control of _2,a^E and _2,b^E follow from the arguments given in the proof of <cit.>. If we sum the obtained bounds over all elements in E in _h and apply a suitable Hölder's inequality, we arrive at

    II_2 ≤𝔅(,T)h^s|E_h|_1,Ø.

A collection of the bounds (<ref>) and (<ref>) shows that

    II≤ℭ C_est/C_data
    |e_h|_1,Ø|E_h|_1,Ø
    +
    𝔅(u,T)h^s|E_h|_1,Ø


Step 2.3. We have now arrived at the estimation of term . The bound for this term follows from the arguments in <cit.>:


    ≲
    𝔊(g) h^s+2|E_h|_1,Ø,


Step 2.4. Finally, if we substitute the bounds for , , and  obtained in (<ref>), (<ref>) and (<ref>), respectively, into (<ref>) we obtain

    |E_h|_1,Ø≤
    C_final,T
    |e_h|_1,Ø
    +
    ℭ(κ,,T)h^s
    +
    𝔊(g) h^s+2.


Step 3. An estimate for e_h. We start with the coercivity bound in (<ref>) for a_h(·;·,·), the definition of e_h, namely e_h = u_I - u_h, and add and subtract the term a(T;,e_h) to obtain

    α_*ν_*|e_h|_1,Ø^2
       ≤ a_h(T_h;e_h,e_h) = a_h(T_h;_I,e_h)-a_h(T_h;_h,e_h)
    
       
    =
    a_h(T_h;_I,e_h)
    -
    a(T;,e_h)
    +
    a(T;,e_h)
    -
    a_h(T_h;_h,e_h).

We now use the first equations of the continuous and discrete systems, (<ref>) and (<ref>), respectively, to obtain

    α_*ν_*|e_h|_1,Ø^2
    ≤
    [
    a_h(T_h;_I,e_h)
    -
    a(T;,e_h)
    ]
    +
    [
    c_N,h^skew(_h;_h,e_h)
    -
    c_N(;,e_h)
    ]
    
    
    +
    [
    c_F,h(_h;_h,e_h)
    -
    c_F(;,e_h)
    ]
    +
    [d_h(_h,e_h)-d(,e_h)]
    +
    (f-f_h,e_h)_0,Ω,

where we have also used that e_h = u_I - u_h ∈_h ⊆𝐙; see Remark <ref>.

Step 3.1. Define ℑ_1 := [
a_h(T_h;_I,e_h)
-
a(T;,e_h)
]. The control of ℑ_1 follows from <cit.>:

    ℑ_1
    ≤C_lip()/C_sol() |E_h|_1,Ø|e_h|_1,Ø
    +
    𝔇(ν,,T)h^s|e_h|_1,Ø.

Here, the constant C_sol() comes from Prop. <ref>: C_sol()|u|_𝐖_q^1(Ω) < 1, and C_lip() denotes a suitable positive constant.

Step 3.2. Define ℑ_2:= c_N,h^skew(_h;_h,e_h) - c_N(;,e_h).
Note that c^skew_N(;,e_h) = c_N(;,e_h) because ∈; see Remark <ref>. So we rewrite the term ℑ_2 as follows:

    ℑ_2 =
     [c_N,h^skew(_h;_h,e_h)
     -
     c_N,h^skew(;,e_h)]
     +
     [c_N,h^skew(;,e_h)
     -
     c^skew_N(;,e_h)].

Since ∈∩𝐇^s+1(Ω) and e_h ∈, we can apply  <cit.> to obtain

    ℑ_2,b
     :=
     c_N,h^skew(;,e_h)
     -
     c^skew_N(;,e_h) ≲𝔈()h^s|e_h|_1,Ø.

Define ℑ_2,a:= c_N,h^skew(_h;_h,e_h)
 -
 c_N,h^skew(;,e_h). Apply <cit.> to obtain

    ℑ_2,a≤C_N C_est/C_data|e_h|^2_1,Ø
     +
     C_N
     (
     C_est/C_data
     +
     C_est/C_data)
     𝔈()h^s|e_h|_1,Ø.

If we combine the previously derived bounds, we arrive at

    ℑ_2
    ≤C_N C_est/C_data
    |e_h|^2_1,Ø
     +
     𝔈()h^s|e_h|_1,Ø.


Step 3.3. Define ℑ_3:= c_F,h(_h;_h,e_h)
-
c_F(;,e_h). We add and subtract c_F,h(;,e_h) to obtain ℑ_3 = [c_F,h(_h;_h,e_h)- c_F,h(;,e_h)] + [c_F,h(;,e_h) - c_F(;,e_h)] =: ℑ_3,a + ℑ_3,b. We bound ℑ_3,a with the help of <cit.>:

    ℑ_3,a≤
    C_for(𝔈()h^s|e_h|_1,Ø + |e_h|^2_1,Ø),

where C_for = C [ (C_est/C_data)^r-2 +  (C_est/C_data)^r-2] and C is the constant in <cit.>. Since u∈∩𝐇^s+1(Ω), a bound for ℑ_3,b follows from <cit.>:

    ℑ_3,b≤
    Ch^s (|u|_s+1,Ω + |u|_s,Ω) |u|^r-2_1,Ø|e_h|_1,Ø≤𝔈()h^s|e_h|_1,Ø.

If we combine the bounds (<ref>) and (<ref>), we obtain

    ℑ_3≤ C_for |e_h|^2_1,Ø + 𝔈()h^s|e_h|_1,Ø.


Step 3.4. Define ℑ_4 := d_h(_h,e_h)-d(,e_h). The control of ℑ_4 is standard. To derive an estimate, we first analyze the local term ℑ_4^E := d_h^E(_h,e_h)-d^E(,e_h). In view of (^0,E_k_h, ^0,E_ke_h - e_h)_0,E = 0, we rewrite ℑ_4^E as follows:

    ℑ_4^E
     =
     (^0,E_k_h - , ^0,E_ke_h)_0,E
     +
     (, ^0,E_ke_h - e_h)_0,E
    
     =
     (^0,E_k_h - , ^0,E_ke_h)_0,E
     +
     ( - ^0,E_k_h, ^0,E_ke_h - e_h)_0,E
     =
     (^0,E_k_h - ,e_h)_0,E.

As a result,
ℑ_4^E
≤^0,E_k_h - _0,Ee_h _0,E.
We control ^0,E_k_h - _0,E as follows:

    ^0,E_k_h - _0,E≤_h - _0,E + ^0,E_k - _0,E
    ≤ - _I _0,E + e_h _0,E
    +
    ^0,E_k - _0,E.

A basic bound for ^0,E_k - _0,E and an application of Lemma <ref> thus show that

    ℑ_4^E
    ≤e_h ^2_0,E
    +
    Ch_E^s+1 |u |_s+1,D(E)e_h _0,E.

If we sum over all the elements in E in _h and apply Hölder's inequality for sums, we obtain an estimate for ℑ_4:

    ℑ_4
    ≤^2 |e_h|_1,Ω^2
     +
    𝔈()h^s+1|e_h|_1,Ω,

where we have also used the Poincaré inequality (<ref>).

Step 3.5. Define ℑ_5 := (f-f_h,e_h)_0,Ω. An estimate for ℑ_5 is direct:

    ℑ_5
    ≲𝔉(f) h^s+2 | e_h |_1,Ω.


Step 3.6. A final estimate for |e_h|_1,Ω. Replace the estimates (<ref>), (<ref>), (<ref>), (<ref>), and (<ref>) into the bounded derived for |e_h|_1,Ø^2 to obtain

    ( α_*ν_* - C_N C_est/C_data - C_for - ^2 )
    |e_h|^2_1,Ø≤C_lip()/C_sol() |E_h|_1,Ø|e_h|_1,Ø
    
    +
    𝔇(ν,,T) h^s|e_h|_1,Ø +
    𝔉(f)h^s+2 |e_h|_1,Ø,

We now replace the bound (<ref>) for |E_h|_1,Ø into the previous estimate to obtain

    |e_h|_1,Ø≤ C_final,( 𝔐(ν,κ,,T) h^s + 𝔑(f,g)h^s+2).


Step 4. A final estimate for |E_h|_1,Ω. We replace (<ref>) into the bound (<ref>) to finally arrive at

    |E_h|_1,Ω≤ C_final,TC_final,u( 𝔐(ν,κ,,T) h^s + 𝔑(f,g)h^s+2)
    
    
    +
    ℭ(κ,,T)h^s
    +
    𝔊(g) h^s+2≤𝔐(ν,κ,,T) h^s + 𝔑(f,g)h^s+2.


Step 5. An estimate for _h. In this last step, we derive an error estimate for the pressure error _h_L^2(Ω). To do so, we let _h ∈_h and use the definition of _h, namely _h = _I - _h and the first equations of the continuous and discrete systems (<ref>) and (<ref>), respectively, to obtain

    b(_h,_h)=b(_h, _I)-b(_h, _h)
    =
    b(_h, _I -)
    +
    b(_h, )-b(_h, _h)
    =
    b(_h, _I- )
    
    
    +
    [a_h(T_h;_h,_h)-a(T;,_h)]
    +
    [c_N,h^skew(_h;_h,_h)-c_N^skew(;,_h)]
    
    
    +
    [c_F,h(_h;_h,_h)-c_F(;,_h)]
    +
    [d_h(_h,_h)-d(,_h)]
    +
    (f-f_h),_h)_0,Ø.


Step 5.1. Define 𝔎_1:=[a_h(T_h;_h,_h)-a(T;,_h)]. The control of the term 𝔎_1 can be found in <cit.>:

    𝔎_1 ≲ |-_h|_1,Ø|_h|_1,Ø
    +
    |T-T_h|_1,Ø|_h|_1,Ø
    +
    𝔇(ν,u,T)h^s|_h|_1,Ø
    ≤𝔐(ν,κ,,T)h^s|_h|_1,Ø
    +
    𝔑(f,g)h^s+2|_h|_1,Ø,

where we have used (<ref>) and (<ref>) to obtain the last bound.

Step 5.2. Define 𝔎_2:= c_N,h^skew(_h;_h,_h)-c_N^skew(;,_h). An estimate for the term 𝔎_2 can be obtained as follows. First, we rewrite 𝔎_2 as

    𝔎_2
    =
    [
    c_N,h^skew(_h;_h,_h)
    -
    c_N,h^skew(;,_h)
    ]
    +
    [
    c_N,h^skew(;,_h)
    -
    c_N^skew(;,_h)]
    
    
    =
    [
    c_N,h^skew(_h - ;_h,_h)
    -
    c_N,h^skew(; - _h,_h)
    ]
    +
    [
    c_N,h^skew(;,_h)
    -
    c_N^skew(;,_h)].

Define
𝔎_2,a:=
c_N,h^skew(_h - ;_h,_h)
-
c_N,h^skew(; - _h,_h). A bound for 𝔎_2,a follows from the use of the bound (<ref>) and the error estimate (<ref>). In fact, we have

    𝔎_2,a   ≤C_N( |_h|_1,Ø + ||_1,Ø )
    | - _h|_1,Ø
    |_h|_1,Ø
       ≤𝔐(ν,κ,,T)h^s|_h|_1,Ø
    +
    𝔑(f,g)h^s+2|_h|_1,Ø.

It thus suffices to bound 𝔎_2,b:= c_N,h^skew(;,_h) - c_N^skew(;,_h). Since ∈𝐇^s+1(Ω) ∩, we can use the bound in <cit.> and obtain
𝔎_2,b≤𝔈()h^s|_h|_1,Ø.
This bound and (<ref>) yield the control of 𝔎_2:

    𝔎_2≤𝔐(ν,κ,,T)h^s|_h|_1,Ø
    +
    𝔑(f,g)h^s+2|_h|_1,Ø.


Step 5.3. Define 𝔎_3:= c_F,h(_h;_h,_h)-c_F(;,_h). We add and subtract the term c_F,h(;,_h) and rewrite 𝔎_3 as follows:

    𝔎_3 = [c_F,h(_h;_h,_h)-c_F,h(;,_h)]
    +
    [c_F,h(;,_h)-c_F(;,_h)].

Define 𝔎_3,a = c_F,h(_h;_h,_h)-c_F,h(;,_h). A direct application of the bound in <cit.> shows that

    𝔎_3,a≤
    C(|_h|_1,Ø^r-2 + ||_1,Ø^r-2)|-_h|_1,Ø|_h|_1,Ø≤
    C_for|-_h|_1,Ø|_h|_1,Ø
    ≤𝔐(ν,κ,,T)h^s|_h|_1,Ø
    +
    𝔑(f,g)h^s+2|_h|_1,Ø.

Define 𝔎_3,b :=c_F,h(;,_h)-c_F(;,_h). Since ∈∩𝐇^s+1(Ω), we can directly apply the bound from <cit.> to obtain

    𝔎_3,b≤ C h^s(||_s,Ø+||_1+s,Ø)||_1,Ø^r-2|_h|_1,Ø≤𝔈() h^s |_h|_1,Ø.

The two estimates derived above allow us to control the term 𝔎_3:

    𝔎_3
    ≤𝔐(ν,κ,,T)h^s|_h|_1,Ø
    +
    𝔑(f,g)h^s+2|_h|_1,Ø.


Step 5.4. Define 𝔎_4 := d_h(_h,_h)-d(,_h) and 𝔎_4^E := d_h^E(_h,_h)-d^E(,_h) for E ∈_h. Following the arguments that lead to (<ref>), we deduce that

    𝔎_4^E =
    (^0,E_k_h - ,v_h)_0,E   ≤^0,E_k_h - _0,E |v_h|_1,E
       ≤
    (
    _h - _0,E
    +
    ^0,E_k  - _0,E
    )
    |v_h|_1,E.

If we sum over all the elements E ∈_h and apply Hölder's inequality, we obtain

    𝔎_4
       ≲
     (
    | _h -  |_1,Ω
    +
    ^0,E_k  - _0,Ω
    )
    |v_h|_1,Ω
       ≤𝔐(ν,κ,,T)h^s|_h|_1,Ø
    +
    𝔑(f,g)h^s+2|_h|_1,Ø.


Step 5.5. Define 𝔎_5 := (f-f_h,_h)_0,Ø and 𝔎_6 := b(_h, _I - ). The term 𝔎_5 was already estimated in (<ref>). A bound for the term 𝔎_6 follows from the definition of _I and the Bramble–Hilbert bound (<ref>): 𝔎_6 ≲𝔓() h^s |_h|_1,Ø.

Step 5.6. In view of the discrete inf-sup condition (<ref>), the bounds (<ref>), (<ref>), (<ref>), and (<ref>), we conclude

    - _h_0,Ø≤𝔒(ν,κ,,T,)h^s
    +
    𝔑(f,g)h^s+2.

This concludes the proof.




§ NUMERICAL EXPERIMENTS



In this section, we report some numerical experiments to evaluate the performance of the proposed VEM. Our goal is to compute the experimental convergence rates in the norms used in the theoretical analysis. The results of this section were obtained using a MATLAB code with k=2, where the nonlinear problem (<ref>) was solved using the fixed-point iteration described in Algorithm 1. The refinement parameter N used to characterize each mesh is the number of elements on each edge of Ω. In Figure <ref> we present examples of the meshes we will use for our tests.



 §.§ Fixed-point iteration

Let us describe the fixed-point iteration used to solve the coupled problem (<ref>). 


Algorithm 1: Fixed-point iteration

13cm0.4pt 

Input: Initial mesh _h, initial guess (_h^0,_h^0,T_h^0)∈_h ×_h×_h, the coefficients ν(·) and κ(·), f_h ∈ [Ł^2(Ø)]^2, g_h∈Ł^2(Ø), r∈[3,4], and =10^-6.


  1: For n≥0, find (_h^n+1,_h^n+1)∈_h×_h such that

    a_h(T_h^n;_h^n+1,_h) + c_N,h^skew(_h^n;_h^n+1,_h)
    +
    c_F,h(_h^n;_h^n+1,_h) 
    
    +
    d_h(_h^n+1,_h)
    +
    b(_h,_h^n+1)   = (f_h,_h)_0,Ø,
    
    
    b(_h^n+1,_h)    = 0,

for every (_h,_h)∈_h×_h. Then, T_h^n+1∈_h is found as the solution of

    𝔞_h(T_h^n;T_h^n+1,S_h)+𝔠_h^skew(_h^n+1;T_h^n+1,S_h)=(g_h,S_h)_0,Ø,

for every S_h∈_h.

  2: If |(_h^n+1,_h^n+1,T_h^n+1)-(_h^n,_h^n,T_h^n)|>, set n← n+1 and go to step 1. Otherwise, return (_h,_h,T_h)=(_h^n+1,_h^n+1,T_h^n+1). Here, |·| denotes the Euclidean norm.

13cm0.4pt



To complete the proposed VEM (<ref>), we need to describe the bilinear forms S_V^E(·,·) and S_T^E(·,·) that satisfy (<ref>) and (<ref>), respectively. For this purpose, we proceed as in <cit.> and use the so-called  dofi–dofi stabilization, which is defined as follows: for each E, we denote by _h, _h and T⃗_h, R⃗_h, the real-valued vectors containg the values of the local degrees of freedom associated with _h, _h in  and with T_h, R_h in _h, respectively. Then, we set

    S_V^E(_h,_h):=_h·_⃗h⃗,
        
    S_T^E(T_h,R_h):=T⃗_⃗h⃗·R⃗_⃗h⃗.


To calculate the error between the exact velocity component  of the solution and the corresponding approximation obtained with our VEM, namely, _h, we will use the following computable quantity

    |e_|_1,Ø:=(∑_E∈_h|-_2^0,E_h|_1,E^2)^1/2.

A similar expression is used to calculate the error between T and T_h. Finally, to calculate the pressure error we consider e__0,Ø:=-_h_0,Ø.



 §.§ Unit square
 Let us first test our method on the unit square Ø=(0,1)^2. The viscosity coefficient ν, the thermal diffusivity coefficient κ, and the parameter r are as follows <cit.>:


  * Test 1: ν(T)=1+T, κ(T)=1+sin(T), and r=3.

  * Test 2: ν(T)=1+e^-T, κ(T)=2+sin(T), and r=4.

The data f and g are chosen so that the exact solution to problem (<ref>) is

    (x_1,x_2)=(-x_1^2(x_1-1)^2x_2(x_2-1)(2x_2-1),x_2^2(x_2-1)^2x_1(x_1-1)(2x_1-1)),

(x_1,x_2)=x_1x_2(1-x_1)(1-x_2)-1/36, and T(x_1,x_2)=x_1^2x_2^2(1-x_1)^2(1-x_2)^2.¨

In Figures <ref> and <ref>, we present the experimental convergence rates in the semi-norm for the velocity error and the temperature error and in the Ł^2-norm for the pressure error for tests 1 and 2. We have calculated these errors for the six families of meshes shown in Figure <ref> with different levels of refinement. The plots also include a reference line with slope -2, which indicates the optimal convergence rate of the method.




Figures <ref> and <ref> show that the theoretical predictions of Theorem <ref> are confirmed. We also note that optimal experimental convergence rates are obtained for each of the mesh families considered. This is computational evidence of the robustness of the VEM with respect to the geometry of the mesh. Figure <ref> shows the solution obtained with our method using a particular family of polygonal meshes.



The next goal is to prove computationally that our VEM is divergence-free. To accomplish this task, we calculate the value of the Ł^2-norm of the discrete divergence of the velocity field _h. In Table <ref>, we present the values of ÷_h_0,Ø that we computed for five mesh families at different refinement levels. As the results in Table <ref> show, the method is indeed divergence-free, and this again is independent of the polygonal mesh used.




 §.§ Influence of the viscosity

In this test, we investigate the influence of the viscosity coefficient ν(·) on the calculation of the error in the approximation of the velocity, pressure, and temperature variables. It is well-known that most standard approximation techniques are not robust with respect to this coefficient. To perform this study, we choose different values of ν and investigate the experimental convergence rates on different polygonal meshes. We use a similar configuration as in the previous tests: Ø=(0,1)^2, s=3, and the solution (, p, T) as previously described. For simplicity, we assume that κ=1 and consider constant values for the viscosity parameter: ν=10^-1, 10^-4, 10^-8.






Our results are shown in Figures <ref>, <ref> and <ref>. In Figure <ref>, we observe that the optimal experimental convergence rate for the velocity error is no longer achieved when the viscosity coefficient becomes smaller. This phenomenon occurs for different meshes, which leads to the conclusion that the loss of convergence for the velocity does not depend on the geometry of the mesh, but only on the physical parameter. A different behavior is observed for the temperature error and the pressure error: the optimal convergence rate is also attained when the viscosity becomes smaller.  We would like to note that this is independent of the considered polygonal mesh. The results obtained for the velocity error and the pressure error are in agreement with the numerical tests reported in <cit.>. Our results also show that the method works for the  temperature error in the same way as for the pressure error. We emphasize that the derivation of an optimal error estimate that is robust with respect to the coefficients ν(·) and κ(·) is not analyzed in our work and is a valuable option for future extensions of the study presented here.


 

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared
to influence the work reported in this paper.

 

Data availability

Data will be made available on request.

siam


