Consistency, convergence and error estimates for a mixed finite element / finite volume scheme to compressible Navier-Stokes equations with general inflow/outflow boundary data
aa r X i v : . [ m a t h . NA ] M a y Consistency, convergence and error estimates for a mixed finiteelement–finite volume scheme to compressible Navier-Stokes equationswith general inflow/outflow boundary data
Young-Sam Kwon ∗ Antonin Novotny † Department of Mathematics, Dong-A UniversityBusan 604-714, Republic of Korea, [email protected] of Toulon, IMATH, EA 2134, BP 20139839 57 La Garde, France, [email protected]
Abstract
We study convergence of a mixed finite element-finite volume scheme for the compressible Navier-Stokes equations in the isentropic regime under the full range 1 < γ < ∞ of the adiabatic coefficient γ for the problem with general non zero inflow-outflow boundary conditions. We propose a modificationof Karper’s scheme [Numer. Math. 125:441-510, 2013] in order to accommodate the non zero boundarydata, prove existence of its solutions, establish the stability and uniform estimates, derive a convenientconsistency formulation of the balance laws and use it to show the weak convergence of the numericalsolutions to a dissipative solution with the Reynolds defect introduced in Abbatiello et al. [Preprint Arxiv1912.12896]. If the target system admits a strong solution then the convergence is strong towards thestrong solution. Moreover, we establish the convergence rate of the strong convergence in terms of thesize of the space discretization h (which is supposed to be comparable with the time step ∆ t ). In the caseof non zero inflow-outflow boundary data, all results are new. The latter result is new also for the no-slipboundary conditions. Key words:
Navier-Stokes equations, Compressible fluids, Dissipative solutions, Defect measure, Crouzeix-Raviart finite element method, Finite volume method, Stability, Convergence, Error estimate
AMS classifiaction:
Contents ∗ The work of the first author was partially supported by NRF- 2017R1D1A1B03030249 and NRF-2019H1D3A2A01101128. † The work of the second author was supported by Brain Pool program funded by the Ministry of Science and ICT throughthe National Research Foundation of Korea, NRF-2019H1D3A2A01101128. Numerical setting 8 L p − L q estimates for finite elements . . . . . . . . . . . . . . . . 164.1.4 Some formulas related to upwinding . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164.2 Preliminaries from mathematical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
11 Concluding remarks 45 Introduction
Evolution of the density r = r ( t, x ) and velocity u = u ( t, x ) through the time interval [0 , T ), T > t ∈ [0 , T ) in a bounded (Lipschitz) domain Ω, x ∈ Ω of a viscous compressible fluid can be described bythe Navier-Stokes equations.
Unlike most of the theoretical literature , in this paper, we consider the Navier-Stokes system with the nonzero inflow-outflow boundary conditions , which is, from the point of view ofapplications, a more appropriate setting, than the ”standard” setting with the no-slip or Navier boudaryconditions. The equations read: ∂ t r + div x ( ru ) = 0 , (1.1) ∂ t ( ru ) + div x ( ru ⊗ u ) + ∇ x p ( r ) = div x S ( ∇ x u ) , S ( ∇ x u ) = µ ∇ x u + ( µ + λ )div u I , µ > , λ + 23 µ > with initial and boundary conditions, r (0) = r , u (0) = u , u | ∂ Ω = u B , r | Γ in = r B , (1.2)where (for the case of simplicity), < r ∈ C (Ω) , u ∈ C (Ω) , (1.3)0 < r B ∈ C (Ω) , u B ∈ C (Ω) (1.4)are given initial and boundary data, andΓ in = ∂ Ω \ Γ out , Γ out = { x ∈ ∂ Ω | u B ( x ) · n ( x ) ≥ n ( x ) is not defined } (1.5)(where n is the outer normal to ∂ Ω), is the inflow boundary. Here, ” { . . . } ” means the closure with respectto the trace topology of R on ∂ Ω.In the above the pressure p is a given function of density r which will be specified later. A typicalexample which will be treated in the paper includes pressure functions p ( r ) ≈ r γ , for large values of r , γ > . In the mathematical literature there is a large variety of numerical methods for solving efficiently com-pressible Navier-Stokes or Euler equations, see [26], [51], [60], [9], [27], [28], [35], [36], [11] – for the finitevolume and/or a combination of finite volume and the finite element methods, e.g. [29], [30] – for the dicon-tinuous Galerkin methods and e.g. [43], [42], [37],[39] - for the finite volume/finite difference schemes – andthe references quoted there. Although these methods give quite convincing results in numerical experimentsand engineering applications, their mathematical analysis is less advanced, see the works of Tadmor et al.[58], [59], [60], Gallouet et al. [35], [36], [12], [38], Jovanovic, Rohde [46], [45] for some pertinent results inthis direction. Namely, the rigorous results about the convergence of the numerical solutions to the (weak, One usually writes S ( ∇ x u ) in its frame indifferent form, S ( ∇ x u ) = µ ( ∇ x u + ( ∇ x u ) T ) + λ div u I . Both writing are equivalentin the strong formulation but not in the weak formulation later and in the numerical scheme The regularity of initial conditions could be relaxed up to0 < r ∈ L γ (Ω) , u ∈ L (Ω) , r u ∈ L (Ω) . γ > γ >
1, see [21], [22], [23], [44], and theproofs, again, are quite scheme structure dependent. Moreover, in these papers, the weak strong uniquenessproperty of the class of Young-measure valued solutions translates to the proofs of strong convergence ofthe numerical solutions to a strong solution, provided the latter exists. This approach, however, does notprovide the rate of convergence in terms of the size of the discretization, the so called error estimates.In 2016, Gallouet et al. [38] established the rigorous error estimates for the Karper’s scheme and theadiabatic exponents γ > /
2. These works were followed by other works establishing error estimates fora few other schemes, see [39], [54], [19], [24]. We remark that the proofs of the error estimates are muchmore ”structure dependent” than the convergence proofs and that they are available only for the adiabaticexponents γ > / γ ≤ / γ > /
2, and it is nowadays a standard result, see [25] (and monographs by Lions [53], Feireisl [13],Novotny, Straskraba [56]) for the no-slip boundary conditions. The same result for the nonzero inflow-outflow boundary conditions is more recent, see [6], [7], [52] preceded by Girinon [40]. This is just to saythat the proof of the convergence of any numerical scheme to a weak solution below the treshold γ = 3 / γ = 3 would be a tremendousprogress. Even for the Karper’s scheme and adiabatic coefficients γ > adapt the Karper’s numerical scheme to the nonzero inflow-outflow boundaryconditions, to prove existence of numerical solutions for this scheme on the polygonal domains and to studythe convergence of its solutions to the dissipative solutions with Reynolds defect introduced in [1] for the whole scale of adiabatic exponents < γ < ∞ . We did not choose Karper’s scheme for its computationalperformance, but rather because it is recognized for its structural closeness to the original system whichallows to adapt to it in a large extend the well known techniques from the ”continuous case”. However,when doing this we have always in mind the universality of the approach and the transportability of the”main steps” of the method to other numerical schemes.4e consider the Karper’s scheme in several versions: 1) The ”classical” version with artificial densitydissipation. 2) Without artificial density dissipation but with artificial pressure perturbation. 3) Combiningboth features mentioned above.Our approach is based on the following steps:1. We modify Karper’s scheme in order to be able to accommodate the non-zero inflow/outflow boundarydata , we prove existence of numerical solutions with the positive density.2. We deduce from the scheme the balance of mass and the balance of energy and we use the balance ofenergy to derive the uniform estimates.3. We rewrite the numerical continuity equation and the numerical momentum equation in their vari-ational forms letting appear ”remainders” vanishing as space h and time ∆ t discretizations tend tozero.4. We rewrite the numerical balance of energy in a consistent form compatible with the variational for-mulations of the balance of mass and momentum .5. Using the consistency formulation of the balance of mass, momentum, and energy , and employingthe uniform estimates, we show that the numerical solution generates a dissipative solution with theReynolds defect in the sense of [1]. This part provides another proof of the existence of the dissipativesolutions via a numerical scheme , and by itself, it is of independent interest.6. Last but not least, by using the weak strong uniqueness principle established in [1] for the limitingdissipative solution we show that the whole sequence of numerical solutions converges strongly to thepointwise classical solution of the original problem provided the latter exists.7. Finally using again the consistency formulations of the continuity and momentum, and employing inaddition the consistency formulation of the energy balance, we derive, by mimicking the continuouscase, the relative energy inequality between the numerical and strong solutions and deduce the errorestimates. All convergence results are ”unconditional” : to obtain them we use only a priory estimates derived fromthe numerical scheme . No a posteriori bounds are needed, in contrast with the most mathematical literatureabout the subject.To the best of our knowledge, all results listed above are new in the case of general non-zero inflow/outflowboundary conditions.
They hold true also in the case of zero boundary conditions , in which situation the errorestimates mentioned in the point 7 is a new result , as well. In the latter case, the convergence to dissipativesolutions with Reynolds stress is reminiscent to the convergence to Young measure valued solutions provedin Feireisl et al. [14]. The advantage of the dissipative solutions with the Reynolds stress may dwell inthe fact that the convergence proofs in this setting require slightly less from the structure of the schemeand thus, may be applicable to a larger class of numerical methods. The pertinence of this statement,however, remains to be confirmed. All results in this paper are obtained through an extensive applicationof the functional analysis and theory of PDEs in the numerical analysis which goes far beyond the standardapproaches current in this domain of research.
The paper is organized as follows. In the Introduction, we fix the structural assumptions, define thedissipative solutions with the Reynolds defect and report the weak-strong uniqueness principle in this class.In Section 2 we describe the numerical setting and introduce the discrete functional spaces related to it.5n Section 3 we suggest the adaptation of Karper’s numerical scheme, see (3.1–3.3), and state the mainresults. Theorem 3.1 about the the existence of numerical solutions and Theorem 3.2 consisting of threeparts: 1) Convergence to dissipative solutions (item 1). 2) Strong convergence to strong solutions (item 2a).3) Error estimates related to the strong convergence (item 2b). Section 4 is reporting necessary numericaland analytical preliminaries for the treatement of the problem. Then, Theorem 3.1 is proved in Section5, and the proof of Theorem 3.2 is established through Sections 6–10. More exactly, we derive energybalance and uniform estimates for the numerical solutions in Section 6 in Lemmas 6.1 and 6.2. In Section7 in Lemmas 7.1, 7.2 and in formulas (7.23–7.24) we derive the convenient consistency formulations ofcontinuity, momentum and energy balance. In Section 8 we use the consistency formulation and uniformestimates to establish the weak convergence to dissipative solutions. This finishes the proof of the first partof Theorem 3.2. In Section 9, we use the weak strong uniqueness principle in the class of dissipative solutionsto prove the strong convergence to a strong solution. This concludes the proofs of part 2 of Theorem 3.2.Finally, in Section 10, we establish the relative energy inequality between the numerical and strong solutionsand conclude the proof of part 3 of Theorem 3.2 by Lemma 10.3 which states the error estimates. We finishthe paper by several concluding remarks in Section 11.We conclude this introductory part by a remark concerning the notation: The special functional spacesare always defined in the text. For the classical Lebesgue, Sobolev, Bochner spaces and their duals, we usethe standard notation, see e.g. Evans [10]. Strong convergence in a Banach space is always denoted ” → ”,while ” ⇀ ” means the weak convergence and ” ⇀ ∗ ” means the star-weak convergence. We suppose p ∈ C [0 , ∞ ) , p (0) = 0 , p ′ ( ̺ ) > ̺ > p by zero to the negative real line so that p ∈ C ( R ) if needed). We associate to p itsHelmholtz function H , H ( ̺ ) = ̺ Z ̺ p ( z ) z d z in particular ̺H ′ ( ̺ ) − H ( ̺ ) = p ( ̺ ) . (1.7)and we shall additionally suppose that H − ap, ap − H are convex functions for some 0 < a < a. (1.8)It is to be noticed that an iconic example of the isentropic pressure p ( ̺ ) = a̺ γ , a > γ > B ∈ C ([0 , ∞ )), we define E B ( ̺ | r ) = B ( ̺ ) − B ′ ( r )( ̺ − r ) − B ( r )and we write E H = E for the particular B = H . For the further use, we also introduce the relative energyfunctional E ( ̺, u | r, U ) = Z Ω (cid:16) ̺ | u − U | d x + E ( ̺ | r ) (cid:17) d x (1.9)to measure a ”distance” between vector fields u , U and positive scalar fields ̺ and r . We remark that (1.7) in combination with (1.6) yields, in particular, H (0) = 0. .3 Dissipative solutions with Reynolds defect Definition 1.1 [Dissipative solution with Reynolds defect]
The quantity [ r , u ] is a dissipative solution of the problem (1.1–1.5) in (0 , T ) × Ω if the following issatisfied:1. ≤ r ∈ C weak ([0 , T ]; L γ (Ω)) ∩ L γ (0 , T ; L γ ( ∂ Ω; | u B · n | d S x )) , for some γ > , v = u − u B ∈ L (0 , T ; W , (Ω; R )) , m ∈ C weak ([0 , T ]; L γγ +1 (Ω; R )) , m = ru a.e. in (0 , T ) × Ω , ru ∈ L ∞ (0 , T ; L (Ω)) , p ( r ) ∈ L ((0 , T ) × Ω)); (1.10)
2. The continuity equation (cid:20)Z Ω r ϕ d x (cid:21) t = τt =0 + Z τ Z Γ out ϕ ru B · n d S x + Z τ Z Γ in ϕ r B u B · n d S x (1.11)= Z τ Z Ω h r ∂ t ϕ + ru · ∇ x ϕ i d x d t, r (0 , · ) = r , holds for any ≤ τ ≤ T , and any test function ϕ ∈ C ([0 , T ] × Ω) ;3. There is a tensor R ∈ L ∞ (0 , T ; M + (Ω; R × )) , where M + (Ω; R × ) denotes the space of positively semi-definite tensor valued (Radon) measures on Ω , such that the integral identity (cid:20)Z Ω m · φ d x (cid:21) t = τt =0 = Z τ Z Ω h ru · ∂ t φ + ru ⊗ u : ∇ x φ + p ( r )div x φ − S ( ∇ x u ) : ∇ x φ i d x d t + Z τ Z Ω ∇ x φ : d R ( t ) d t, m (0 , · ) = r u (1.12) holds for any ≤ τ ≤ T and any φ ∈ C ([0 , T ] × Ω; R ) , φ | ∂ Ω = 0 ;4. The energy inequality Z Ω (cid:20) r | v | + H ( r ) (cid:21) ( τ ) d x + D Z Ω dTr[ R ( τ )] + Z τ Z Ω S ( ∇ x u ) : ∇ x u d x d t (1.13)+ Z τ Z Γ out H ( r ) u B · n d S x d t + Z τ Z Γ in H ( r B ) u B · n d S x d t ≤ Z Ω (cid:20) r | v | + H ( r ) (cid:21) d x − Z τ Z Ω [ ru ⊗ u + p ( r ) I ] : ∇ x u B d x d t + Z τ Z Ω ru · ∇ x u B · u B d x d t + Z τ Z Ω S ( ∇ x u ) : ∇ x u B d x d t − Z τ Z Ω ∇ x u B : d R ( t )d t holds for a.a. τ ∈ (0 , T ) with D = min { / , a/ } , v = u − u B . We say that r ∈ C weak ([0 , T ]; X ), X a Banach space, if r : [0 , T ] X is defined on [0 , T ], r ∈ L ∞ (0 , T ; X ) and < F , r > X ∗ ,X ∈ C [0 , T ] with any F in the dual space X ∗ to X . This means that for any ξ ∈ R and a. a. t ∈ (0 , T ), ξ T R ( t ) ξ is a positive Radon measure on Ω and ∀ ϕ ∈ C (Ω), R Ω ϕ ( x )d R ( · ) ∈ L ∞ ((0 , T )). Lemma 1.1 [Weak–strong uniqueness for dissipative solutions]
Let Ω ⊂ R be a bounded Lipschitz domain.Suppose that p satisfies the hypotheses (1.6–1.8) and the initial and boundary conditions verify (1.3–1.4).Let [ r , u ] be a dissipative solution of problem (1.1–1.2) in the sense of Definition 1.1, and let [ r, U = V + u B ] be a strong solution of the same problem belonging to the class V ∈ C ([0 , T ] × Ω; R ) , ∇ x V ∈ C ([0 , T ] × Ω; R ) , V | ∂ Ω = 0 , r ∈ C ([0 , T ] × Ω) , r := inf (0 ,T ) × Ω r > emanating from the same initial and boundary data. Then r = r, u = U in (0 , T ) × Ω , R = 0 . We suppose that the physical space is a polyhedral bounded domain Ω ⊂ R that admits a tetrahedral mesh T = T h ; the individual elements in the mesh will be denoted by K = K h ∈ T (closed sets) and their gravitycenters by x K . Faces in the mesh are denoted as σ = σ h (close sets in R ) and their gravity centers by x σ ,whereas E = E h is the set of all faces. We also denote by E ( K ) the set of all faces of K ∈ T . The set of faces σ ⊂ ∂ Ω is denoted E ext , while E int = E \ E ext . Finally, we suppose that the mesh fits to the inflow-outflowboundaries, meaning that E ext = E in ∪ E out , Γ in = ∪ σ ∈E in σ, Γ out = ∪ σ ∈E out σ, (2.1)where E in , E out are defined in (2.7) later.We denote by h K the diameter of K and by h K the radius of the largest ball included in K . We call h = sup K ∈T h K the size of the mesh and denote h = inf K ∈T h K .For two numerical quantities a , b , we shall write a < ∼ b if a ≤ cb, c > , a ≈ b if a < ∼ b and b < ∼ a. Here, “constant” typically means a generic quantity independent of the size h of the mesh and the time step∆ t used in the numerical scheme as well as other parameters as the case may be.In addition, we require the mesh to be admissible in the sense of Eymard et al. [11, Definition 2.1]:1. For K, L ∈ T , K = L , the intersection K ∩ L , if non-empty, is either a vertex, or an edge, or a face σ ∈ E . In the latter case, we write σ = K | L .2. There holds h ≈ h We denote by n σ,K the unit normal to the face σ ∈ E ( K ) outwards to K . On the other hand we associateto each element σ ∈ E a fixed normal vector n = n σ . If σ ∈ E ext then n σ is always the outer normal to ∂ Ω. In the sequel we shall omit in the notation the dependence on the ”size” h whenever there is no danger of confusion. .2 Piecewise constant finite elements We introduce the space Q (Ω) = Q h (Ω) = n g ∈ L (Ω) (cid:12)(cid:12)(cid:12) g | K = a K ∈ R o of piecewise constant functions along with the associated projectionΠ Q = Π Qh : L (Ω) → Q (Ω) , b g | K := Π Q [ g ] | K = g K := 1 | K | Z K g d x. For a function g ∈ Q (Ω) and any σ ∈ E int , we denote g + σ := g + n σ = lim δ → g ( x σ + δ n ) , g − σ := g n σ = lim δ → g ( x σ − δ n ) . (2.2)Further, we define the jumps and mean values over σ (relative to n σ ),[[ g ]] σ = [[ g ]] σ, n := g + σ − g − σ , { g } σ := 12 (cid:0) g + + g − (cid:1) . (2.3) A differential operator D acting on the x − variable will be discretized as D h v | K = D ( v | K ) for any v differentiable on each element K ∈ T . The
Crouzeix-Raviart finite element spaces (see Brezzi and Fortin [5], among others) are defined as V (Ω) = V h (Ω) := n v ∈ L (Ω) (cid:12)(cid:12)(cid:12) v | K = affine function , (2.4) Z σ v | K dS x = Z σ v | L dS x for any σ = K | L ∈ E int o , together with V (Ω) = V h, = (cid:26) v ∈ V (Ω) (cid:12)(cid:12)(cid:12) Z σ v dS x = 0 for any σ ∈ E ext (cid:27) . (2.5)Next, we introduce the associated projectionΠ V = Π Vh : W , (Ω) → V (Ω) , ˜ v = Π V [ v ] := X σ ∈E v σ φ σ , (2.6)where v σ = 1 | σ | Z σ v d S x and { φ σ } σ ∈E ⊂ V (Ω) is a basis in V (Ω) given by1 | σ ′ | Z σ ′ φ σ = δ σ,σ ′ , ( σ, σ ′ ) ∈ E . .4 Convective terms, upwinds Suppose, that we are given u B ∈ V (Ω , R ). We define E in = { σ ∈ E ext | u B,σ · n < } , E out = { σ ∈ E ext | u B,σ · n ≥ } (2.7)We define for any σ ∈ E int for any g ∈ Q (Ω) its upwind value on any σ ∈ E int , g up σ = (cid:26) g K if σ = K | L ∈ E int and u σ · n σ,K ≥ ,g L if σ = K | L ∈ E int and u σ · n σ,K < (cid:27) . (2.8)Finally, we associate to any face σ ∈ E int the upwind operator Up σ [ g, u ] := Up σ, n [ g, u ] defined asUp σ, n [ g, u ] = g − [ u σ · n ] + + g + [ u σ · n ] − , where [ c ] + = max { c, } , [ c ] − = min { c, } , (2.9)and to any face σ ∈ E ( K ) the specific flux F σ,K (outwards the element K ) defined as F σ,K [ g, u ] = g up σ u σ · n σ,K . (2.10) For simplicity, we shall consider the constant time step ∆ t > T = N ∆ t , N ∈ N and we set I k = ( τ k − , τ k ] , τ k = k ∆ t, k ∈ Z . (2.11)Suppose that we have functions v k : Ω → R , k = 0 , . . . , N . For convenience, we set v k ( x ) = v ( x ) if k ≤ v k ( x ) = v N ( x ) if k > N , and we introduce numbers D t v k ( x ) = v k ( x ) − v k − ( x )∆ t , k ∈ Z . Finally, we define v ( t, x ) = X k ∈ Z I k ( t ) v k ( x ) , D t v ( t, x ) = X k ∈ Z I k ( t ) D t v k ( x ) , (2.12) v e ( t, x ) = X k ∈ Z I k ( t ) (cid:16) v k − ( x ) + ( t − ( k − t ) D t v k ( x ) (cid:17) , so that ∂ t v e ( t, x ) = D t v ( x ). (2.13)In the sequel, we denote by L ∆ t (0 , T ; Q h (Ω)) := L (0 , T ; Q (Ω)) resp. L ∆ t (0 , T ; V h (Ω)) := L (0 , T ; V (Ω)) thespaces of piecewise constant functions from [0 , T ] to Q (Ω) and V (Ω), respectively (constant on each I k andextended by the value in I to the negative real axes and by the value in I N to [ T, ∞ )). Having collected the necessary preliminary material, we are ready to introduce the numerical scheme tosolve the compressible Navier-Stokes system with the nonhomogenous boundary conditions.10 .1 Numerical scheme
We are given the approximations of the initial and boundary conditions ̺ = ̺ h = Π Qh [ r ] , u = u h = Π V [ u ] , ̺ B = ̺ B,h = Π Q [ r B ] , u B = u B,h = Π V [ u B ] (3.1)We are searching for ̺ h, ∆ t ( t, x ) = N X k =1 X K ∈T I k ( t )1 K ( x ) ̺ kK,h, ∆ t , u h, ∆ t ( t, x ) = N X k =1 X K ∈T I k ( t )1 K ( x ) u kK,h, ∆ t where ̺ k ∈ Q (Ω) , ̺ kh, ∆ t > , u kh, ∆ t ∈ V (Ω; R ) , v k = u k − u B ∈ V (Ω; R ) , k = 1 , . . . , N (3.2)such that the following algebraic equations (for the unknowns ̺ kK , u kσ , k = 1 , . . . , N , K ∈ T , σ ∈ E ) aresatisfied: Approximation of the continuity equation Z Ω D t ̺ k φ d x + X K ∈T X σ ∈E ( K ) ∩E int Z σ F σ,K ( ̺ k , u k ) φ d S x + X σ ∈E out Z σ ̺ k u B,σ · n σ φ d S x (3.3)+ X σ ∈E in Z σ ̺ B u B,σ · n σ φ d S x + κh ω X σ ∈E int Z σ [[ ̺ k ]] σ [[ φ ]] σ d S x = 0for all φ ∈ Q (Ω), where κ = 1 or κ = 0 and ω > Approximation of the momentum equation Z Ω D t ( ̺ k b v k ) · b φ d x + X K ∈T X σ ∈E ( K ) ∩E int Z σ F σ,K ( ̺ k b v k , u k ) · b φ d S x (3.4)+ X K ∈T X σ ∈E ( K ) Z σ ̺ k b u k · n σ,K u B · b φ d S x + X σ ∈E out Z σ ̺ k u B,σ · n σ b v k · b φ d S x + X σ ∈E in Z σ ̺ B u B,σ · n σ b v k · b φ d S x + Z Ω (cid:16) S ( ∇ h u k ) : ∇ h φ − p h ( ̺ k )div h φ (cid:17) d x + κh ω X σ ∈E int Z σ [[ ̺ k ]] σ { b v k } σ [[ b φ ]] σ d S x = 0for any φ ∈ V (Ω; R ), where p h ( ̺ ) = p ( ̺ ) + ˜ κh η ̺ , ˜ κ ∈ { , } , η > . (3.5) Remark 3.1
1. It is to be noticed that the background linear momentum ̺ u B in the momentum equationis not ”upwinded”. If it were ”upwinded” we would loose the derivation of uniform estimates from theenergy balance. This observation seems to have an universal character valid for the discretizations ofthe inflow/outflow problems via the finite volume methods, in general. In what follows, we omit the indexes “ h ” and/or “∆ t ” and write simply ̺ k instead of ̺ kh, ∆ t , ̺ instead of ̺ h, ∆ t , etc., in orderto avoid the cumbersome notation, whenever there is no danger of confusion. In agreement with (1.1), here and in the sequel, R Ω S ( ∇ h u k ) : ∇ h φ d x means exactly R Ω ( µ ∇ h u k : ∇ h φ + ( µ + λ )div h u k div h φ ) d x. This form is important for the estimates, since the Korn inequality does not hold in the Crouzeix-Raviartfinite element space. . The term κh ω P σ ∈E int R σ [[ ̺ k ]] σ [[ φ ]] σ d S x and the last term at the last line of the momentum equationcreate an artificial diffusion in the density field. The perturbation ˜ κh η ̺ represents a ”regularization”of pressure. Under certain assumptions on η , the pressure regularization itself guarantees alreadythe convergence of the scheme to the dissipative solutions with Reynolds defect (with or without theartificial density diffusion) for the whole range < γ < ∞ . We have chosen this particular setting tostate the convergence theorem (Theorem 3.2). However, under certain hypotheses on p and for someranges of γ ’s, the artificial pressure perturbation and even the artificial diffusion are not needed. Weshall discuss these situations in the Remark 6.2 after the main theorem. The first result postulates the existence of numerical solutions.
Theorem 3.1 [ Existence of numerical solutions.]
Under assumptions (1.6–1.8), (1.3–1.4), the algebraicsystem (3.1–3.4) admits at least one solution [ ̺ h, ∆ t , u h, ∆ t ] ∈ L ∆ t (0 , T ; Q h (Ω)) × L ∆ t (0 , T ; V h (Ω)) . Any ofits solution has a strictly positive density. Theorem 3.1 will be proved in Section 5. Let us notice that uniqueness of numerical solutions to thisscheme is not known.The main result deals with the case h ≈ ∆ t . It guarantees a (weak) convergence of a subsequence ofnumerical solutions to a dissipative solution with Reynolds defect. Moreover, if the original problem admitsa strong solution, than the weak convergence of a subsequence becomes a strong convergence of the wholesequence to the strong solution. Even more, in the latter case the convergence rate of the strong convergence can be evaluated in terms of a power of the discretization parameter h : in the other words we have an errorestimate. Theorem 3.2 [Main theorem: Convergence and error estimates.]
Let h = ∆ t and ˜ κ = 1 , and η ∈ (0 , / if κ = 0 or η ∈ (0 , min { ω, / } ) if κ = 1 . (3.6) Suppose that the pressure satisfies assumptions (1.6–1.8) and that the initial and boundary conditions verify(1.3–1.4). Consider a sequence of numerical solutions [ ̺ h , u h ] ∈ L h (0 , T ; Q h (Ω)) × L h (0 , T ; V h (Ω)) of theproblem (3.1–3.4). Then we have:1. There exists a subsequence [ ̺ h , u h ] (not relabeled) such that ̺ h ⇀ r in L ∞ (0 , T ; L γ (Ω)) , γ = 1 + 1 /a, (3.7) v h ⇀ v in L (0 , T ; L (Ω)) , (3.8) ̺ h b u h ⇀ ∗ ru in L ∞ (0 , T ; L γγ +1 (Ω)) (3.9) (where ru ⊗ u , p ( r ) ∈ L ∞ (0 , T ; L (Ω)) , u = v + u B ) and ̺ h b u h ⊗ b u h + p ( ̺ h ) I ⇀ ∗ ( ru ⊗ u + p ( r ) I ) in L ∞ (0 , T ; M (Ω; R × )) (3.10) where [ r , u ] and R = ( ru ⊗ u + p ( r ) I ) − ( ru ⊗ u + p ( r ) I ) is a dissipative solution of the problem (1.1–1.5) in the sense of Definition 1.1. To simplify notation, we denote ̺ h,h by ̺ h or even ̺ , u h,h by u h or even u etc., if there is no danger of confusion. . If the problem (1.1–1.5) admits a strong solution [ r, U ] in the class (1.14) then there holds:(a) The (weak) limit [ r , u ] in (3.7–3.8) is equal to [ r, U = V + u B ] and the defect R = 0 . In this casethe whole sequence [ ̺ h , u h ] converges to [ r, U ] in the sense (3.7–3.9), and, moreover, p h ( ̺ h ) ⇀ p ( r ) in L ((0 , T ) × Ω) , (3.11) ̺ h b u h ⇀ r U in L ((0 , T ) × Ω) . (3.12) (b) There exists α = α ( η ) > (in the case κ = 0 ) and α = α ( η, ω ) > (in the case κ = 1 ), and apositive number c dependent on r := inf (0 ,T ) × Ω r, r := sup (0 ,T ) × Ω r, k∇ x r, ∂ t r, V , ∇ x V , ∇ x V , ∂ t V , u B , ∇ x u B , ∇ x u B k C ([0 ,T ] × Ω) (3.13) such that h E (cid:16) ̺ h , b v h | r, V (cid:17)i τ + Z τ (cid:16) k u h − ˜ U k L (Ω) + k∇ h ( u h − ˜ U ) k L (Ω) (cid:17) d t < ∼ c h α . (3.14) Remark 3.2
1. Theorem 3.2 holds, in particular, for the isentropic pressure p ( ̺ ) = a̺ γ , a > , < γ < ∞ .2. Formulas (3.11) and (3.12) imply strong convergence in many particular pertinent situations: Forexample, for the ”iconic” isentropic pressure p ( ̺ ) = a̺ γ , a > we have ̺ h → r in L γ ((0 , T ) × Ω) , √ ̺ h b u h → √ r U in L ((0 , T ) × Ω) for the whole sequence of numerical approximations. Moreover, formula (3.14) implies, in addition, u h → U , ∇ h u h → ∇ x U in L ((0 , T ) × Ω) for the whole sequence.3. Suppose now that ˜ κ = 0 . In this situation, we shall suppose for the pressure p in addition to (1.6–1.8)also ∀ ̺ > , π̺ a − + p̺ γ − ≤ p ′ ( ̺ ) , < a ≤ min { γ, } , γ > , where ≤ π , < p . (3.15) In this setting, Theorem 3.2 continues still to hold in the following situations:(a) If in (3.15), π > , and κ ∈ { , } , ω > , ≤ γ < ∞ . (3.16) In this case, if κ = 0 , the number α > in (3.14) depends solely on γ .(b) If / < γ < , and in addition to (1.6–1.8), (3.15) ̺ p ′ ( ̺ ) ̺ is nonincreasing on (0 , ∞ ) , (3.17) provided κ ∈ { , } , ω > . (3.18) In this case, in (3.14), α = α ( γ ) if κ = 0 and α = α ( γ, ω ) if κ = 1 .We notice that this case includes the isentropic pressure p ( ̺ ) = a̺ γ , a > . c) If / < γ ≤ / and if, in addition to (1.6–1.8), the condition (3.17) holds, provided κ = 1 , ω ∈ (0 , γ − . (3.19) In this case, in (3.14), α = α ( γ, ω ) .We notice that this case includes the isentropic pressure p ( ̺ ) = a̺ γ , a > .4. The error estimates derived in Theorem 3.2 require ∆ t ≈ h , and they may not be optimal. Theycertainly are not optimal in the range γ > / when, at least in the case of no-slip boundary condi-tions, ”better” estimates (which do not require the condition ∆ t ≈ h ) can be derived by using anotherapproach, see Gallouet at al. [38]. For the general non zero inflow/outflow setting, this is still aninteresting open problem. We refer the reader to the paper [54] for the discussion on the ”optimality”of the error estimates.5. Local in time existence of strong solutions to problem (1.1–1.5) notably with non zero inflow outflowboundary conditions is discussed in Valli, Zajaczkowski [61]. Theorem 3.2 will be proved through Sections 6–10. Some hints how to modify the proof of its variantsmentioned in the third item of Remark 3.2 are available in Remarks 6.2, 7.1, 7.2.
In this part, we recall several classical inequalities related to the discrete functional spaces which will beused throughout the paper.
We recall Jensen’s inequalities k b v k L q ( K ) < ∼ k v k L q ( K ) for all v ∈ L q ( K ), 1 ≤ q ≤ ∞ , k v σ k L q ( σ ) < ∼ k v k L q ( σ ) for all v ∈ L q ( σ ), 1 ≤ q ≤ ∞ . (4.1)together with the error estimate ∀ v ∈ W ,q (Ω) , ( k v − b v k L q ( K ) < ∼ h k∇ x v k L q ( K ; R ) ∀ v ∈ W ,q ( K ) k v − b v k L q (Ω) < ∼ h k∇ x v k L q (Ω; R ) ∀ v ∈ W ,q (Ω) ) , ≤ q ≤ ∞ . (4.2)We also recall the Poincar´e type inequalities on the mesh elements, k v − v σ k L q ( K ) < ∼ h k∇ x v k L q ( K ; R ) , ∀ σ ∈ E ( K ) , v ∈ W ,q ( K ) , ≤ q ≤ ∞ , k v − v K k L q ( K ) < ∼ h k∇ x v k L q ( K ; R ) , ∀ K ∈ T , v ∈ W ,q ( K ) , ≤ q ≤ ∞ . (4.3)14 .1.2 Properties of piecewise constant and Crouzeix-Raviart finite elements We report the estimates of jumps on mesh elements, k [[ b v ]] σ = K | L k L q ( σ ) < ∼ h k∇ h v k L q ( K ∪ L ; R ) , ∀ v ∈ V (Ω) , ≤ q ≤ ∞ , k v | K − v | L k L q ( K ∪ L ) < ∼ h k∇ h v k L q ( K ∪ L ; R ) , ∀ σ = K | L, v ∈ V (Ω) , ≤ q ≤ ∞ , (4.4)see Gallouet et al. [36, Lemma 2.2].Further we recall a global version of the Poincar´e inequality on V (Ω) k v − b v k L q (Ω) < ∼ h k∇ h v k L q (Ω) for all v ∈ V (Ω), 1 ≤ q < ∞ , (4.5)along with the global error estimate k v − ˜ v k L q (Ω) + h k∇ h ( v − ˜ v ) k L q (Ω; R ) < ∼ h m k v k W m,q (Ω) (4.6)for any v ∈ W m,q (Ω), m = 1 , , ≤ q ≤ ∞ , see Karper [50, Lemma 2.7] or Crouzeix and Raviart [8].Next, we shall deal with the Sobolev properties of piecewise constant and Crouzeix-Raviart finite ele-ments. To this end we introduce a discrete (so called broken) Sobolev H ,p -(semi)norm on Q (Ω), k g k pQ ,p (Ω) = X σ ∈E int Z σ [[ g ]] pσ h p − dS x , ≤ p < ∞ . (4.7)and a discrete (so called broken) Sobolev H ,p -(semi)norm on V (Ω), k v k pV ,p (Ω) = Z Ω |∇ h v | p d x. (4.8)Related to the Q ,p -norm, we report the following discrete Sobolev and Poincar´e type inequalities: Wehave ∀ g ∈ Q (Ω) , k g k L p ∗ (Ω) < ∼ k g k Q ,p (Ω) + k g k L p (Ω) , (4.9)where 1 ≤ p ∗ ≤ p − p if 1 ≤ p <
3, and 1 ≤ p ∗ < ∞ , if p ≥
3, see Bessemoulin-Chatard et al. [4, Theorem 3].Likewise, we have the Discrete Sobolev inequality, ∀ v ∈ V (Ω) , ( k v k L p ∗ ( K ) < ∼ k∇ x v k L p ( K ) + k v k L p ( K ) , k v k L p ∗ (Ω) < ∼ k v k V ,p (Ω) + k v k L p (Ω) ) , ≤ p ≤ ∞ , (4.10) ∀ v ∈ V (Ω) , k v k L p ∗ (Ω) < ∼ k v k V ,p (Ω) , (4.11)see e.g. [38, Lemma 9.3].Finally, we report the well-known identities for the Crouzeix-Raviart finite elements, Z Ω div h ˜ u w d x = Z Ω div x u w d x for any w ∈ Q (Ω) and u ∈ V (Ω; R ) , (4.12) Z Ω ∇ h v ⊗ ∇ h ˜ ϕ d x = Z Ω ∇ h v ⊗ ∇ x ϕ d x for all v ∈ V h (Ω) , ϕ ∈ W , (Ω) , (4.13)see [50, Lemma 2.11]. 15 .1.3 Trace and “negative” L p − L q estimates for finite elements We start by the classical trace estimate, k v k L q ( ∂K ) < ∼ h /q (cid:0) k v k L q ( K ) + h k∇ x v k L q ( K ; R ) (cid:1) , ≤ q ≤ ∞ for any v ∈ C ( K ) , (4.14)The following can be easily obtained from the previous one by means of the scaling arguments: k w k L q ( ∂K ) < ∼ h /q k w k L q ( K ) for any 1 ≤ q ≤∞ , w ∈ P m , (4.15)where P m denotes the space of polynomials of order m .In a similar way, from the local estimate k w k L p ( K ) < ∼ h (cid:16) p − q (cid:17) k w k L q ( K ) ≤ q < p ≤ ∞ , w ∈ P m , (4.16)we deduce the global version k w k L p (Ω) ≤ ch (cid:16) p − q (cid:17) k w k L q (Ω) ≤ q < p ≤ ∞ , for any w | K ∈ P m ( K ) , K ∈ T . (4.17)In particular, for a piecewise constant function in (0 , T = N ∆ t ], a ( t ) = P Nk =1 a k I k ( t ), where I k = (( k − t, k ∆ t ] we have k w k L p ( I k ) < ∼ (∆ t ) (cid:16) p − q (cid:17) k w k L q ( I k ) ≤ q < p ≤ ∞ , (4.18)and k w k L p (0 ,T ) < ∼ (∆ t ) (cid:16) p − q (cid:17) k w k L q (0 ,T ) ≤ q < p ≤ ∞ . (4.19) The following formulas can be easily verified by a direct calculation, see e.g. [17, Section 2.4]1. Local conservation of fluxes: ∀ σ = K | L ∈ E int , F σ,K [ g, u ] = − F σ,L [ g, u ] , ∀ σ ∈ E int , Up σ, n [ g, u ] = − Up σ, − n [ g, u ] and [[ g ]] n σ = − [[ g ]] − n σ . (4.20)2. For all r, g ∈ Q (Ω), u ∈ V (Ω , R ), X K ∈T r K X σ ∈E int | σ | F σ,K [ g, u ] = − X σ ∈E int Up σ [ g, u ][[ r ]] σ . (4.21)3. For all r, g ∈ Q (Ω), u , u B ∈ V (Ω; R ), u − u B ∈ V (Ω; R ), φ ∈ C (Ω), Z Ω g u · ∇ x φ d x = − X K ∈T X σ ∈E ( K ) ∩E int F σ,K [ g, u ] r (4.22)+ X K ∈T X σ ∈E ( K ) ∩E int Z σ ( r − φ ) [[ g ]] σ, n σ,K [ u σ · n σ,K ] − d S x + X K ∈T X σ ∈E ( K ) Z σ g ( u − u σ ) · n σ,K ( φ − r )d S x + Z Ω ( r − φ ) g div h u d x + X σ ∈E out Z σ g u B,σ · n σ ( φ − r )d S x + X σ ∈E in Z σ g u B,σ · n σ ( φ − r )d S x . .2 Preliminaries from mathematical analysis Here we recall some elements of mathematical analysis needed in the proofs. First is the Shaeffer’s fixedpoint theorem, see e.g. Evans [10, Theorem 9.2.4].
Lemma 4.1 [Shaeffer’s fixed point theorem]
Let
Θ : Z Z be a continuous mapping defined on a finitedimensional normed vector space Z . Suppose that the set { z ∈ Z | z = ΛΘ( z ) , Λ ∈ [0 , } is non empty and bounded. Then there exists z ∈ Z such that z = Θ( z ) . The next result is a version of Div-Curl lemma, see Lemma 8.1 in [1].
Lemma 4.2
Let O = (0 , T ) × Ω , where Ω ⊂ R d , d ≥ is a bounded domain. Suppose that r n ⇀ r weakly in L p ( O ) , v n ⇀ v in L q ( O ) , p > , q > , and r n v n ⇀ w in L r ( O ) , r > . In addition, let ∂ t r n = div x g n + h n in D ′ ( O ) , k g n k L s ( O ; R d ) < ∼ , s > , h n precompact in W − ,z ( O ) , z > , and k∇ x v n k M ( O ; R d ) < ∼ uniformly for n → ∞ . Then w = rv a.a. in O. Finally, the last two lemmas are well known results from convex analysis, see e.g. Lemma 2.11 andCorollary 2.2 in Feireisl [13].
Lemma 4.3
Let O ⊂ R d , d ≥ , be a measurable set and { v n } ∞ n =1 a sequence of functions in L ( O ; R M ) such that v n ⇀ v in L ( O ; R M ) . Let
Φ : R M → ( −∞ , ∞ ] be a lower semi-continuous convex function.Then Φ( v ) : O R is integrable and Z O Φ( v )d x ≤ lim inf n →∞ Z O Φ( v n )d x. Lemma 4.4
Let O ⊂ R d , d ≥ be a measurable set and { v n } ∞ n =1 a sequence of functions in L ( O ; R M ) such that v n ⇀ v in L ( O ; R M ) . et Φ : R M → ( −∞ , ∞ ] be a lower semi-continuous convex function such that Φ( v n ) ∈ L ( O ) for any n ,and Φ( v n ) ⇀ Φ( v ) in L ( O ) . Then Φ( v ) ≤ Φ( v ) a.e. in O. (4.23) If, moreover, Φ is strictly convex on an open convex set U ⊂ R M , and Φ( v ) = Φ( v ) a.e. on O, then v n ( y ) → v ( y ) for a.a. y ∈ { y ∈ O | v ( y ) ∈ U } (4.24) extracting a subsequence as the case may be. Using in the discrete continuity equation test function φ ≈ φB ′ ( ̺ k ), φ ∈ Q (Ω) we obtain after an elementarybut laborious calculation the discrete renormalized continuity equation. This result is formulated in thefollowing lemma. Lemma 5.1 [Renormalized discrete continuity equation]
Let B ∈ C ( R ) . Then Z Ω D t B ( ̺ k ) φ d x + X K ∈T X σ ∈E ( K ) ∩E int Z σ F σ,K ( B ( ̺ k ) , u k ) φ d S x (5.1)+ κh ω X σ ∈E int Z σ [[ ̺ k ]] σ [[ φB ′ ( ̺ k )]] σ d S x − Z Ω φ (cid:16) B ( ̺ k ) − ̺ k B ′ ( ̺ k ) (cid:17) div h u k d x + 1∆ t Z Ω E B ( ̺ k − | ̺ k ) φ d x + X σ ∈E int Z σ E B ( ̺ k, + σ | ̺ k, − σ ) | u kσ · n | φ + dS x + X σ ∈E out Z σ B ( ̺ k ) u B,σ · n σ · φ d S x + X σ ∈E in Z σ E B ( ̺ B | ̺ k ) | u B,σ · n σ | φ d S x = X σ ∈E in Z σ B ( ̺ B ) | u B,σ · n σ | φ d S x for all φ ∈ Q (Ω) . It is enough to show that ̺ k > ̺ k − >
0. (Recall that ̺ B > φ = 1 and B the non negative convex function B ( ̺ ) = max {− ̺, } (resp. any of itsnon negative C ( R ) convex approximations B ε such that B ε ( s ) → B ( s ), B ′ ε → B ′ ( s ) for all s = 0). We get, Z Ω B ε ( ̺ k ) d x ≤ Z Ω B ε ( ̺ k − ) d x + ∆ t (cid:16) X σ ∈E in Z σ B ε ( ̺ B ) | u B,σ · n σ | d S x + Z Ω ( B ε ( ̺ k ) − ̺ k B ′ ε ( ̺ k ))div h u k d x (cid:17) . B ( s ) − sB ′ ( s ) = 0 for all s = 0, B ( ̺ k − ) = 0, B ( ̺ B ) = 0, we deduce from the above that B ( ̺ k ) = 0, i.e. ̺ k ≥ := { x ∈ Ω | ̺ k ( x ) = 0 } = ∪ K ∈ ˜ T K , where ∅ 6 = ˜ T ⊂ T . Taking in the continuity equation(3.3) test function φ = 1 Ω , we get − Z Ω ̺ k − d x = − X K ∈ ˜ T X σ ∈ ∂ Ω ∩E int ∩E ( K ) Z σ ̺ k, + n σ,K [ u σ · n σ,K ] − d S x + X σ ∈ ∂ Ω ∩E in Z σ ̺ B | u B,σ · n σ | d S x , where the left hand side is strictly negative while the right hand side is positive. Consequently Ω = ∅ . Our goal is to show the following: Given 0 < ̺ k − ∈ Q (Ω), v k − ∈ V (Ω; R ) (and the boundary data ̺ B , u B ), the algebraic system (3.3-3.4) admits a solution 0 < ̺ k ∈ Q (Ω) and v k ∈ V (Ω; R ).We first prove the solvability of the linear equation (3.3) with u k = v k + u B , v k ∈ V (Ω; R ) given: Itreads Z Ω ̺ k φ K d x + ∆ t X K ∈T X σ ∈E ( K ) ∩E int Z σ F σ,K ( ̺ k , u k ) φ K d S x +∆ t X σ ∈E out Z σ ̺ k u B,σ · n σ φ K d S x +∆ tκh ω X σ ∈E int Z σ [[ ̺ k ]] σ [[ φ K ]] σ d S x = F K , where { φ K } K ∈T is a basis in Q (Ω) and F K = Z Ω ̺ k − φ K d x − ∆ t X σ ∈E in Z σ ̺ B u B,σ · n σ φ K d S x . It admits a unique solutions ̺ k ∈ Q (Ω) – we denote it by η ( v k )– since the corresponding homogenous system Z Ω ̺ k φ K d x + ∆ t X K ∈T X σ ∈E ( K ) ∩E int Z σ F σ,K ( ̺ k , u k ) φ K d S x +∆ t X σ ∈E out Z σ ̺ k u B,σ · n σ φ K d S x +∆ tκh ω X σ ∈E int Z σ [[ ̺ k ]] σ [[ φ K ]] σ d S x = 0admits a unique solution ̺ k = 0 as we easily show by applying to it the reasoning of Section 5.2 with theconvex function B ( s ) = | s | . We already proved in Section 5.2 that η ( v k ) > R Ω η ( v k ) d x is bounded and that the map η is continuous from V (Ω; R ) to Q (Ω).We now define a map Θ : V (Ω; R ) → V (Ω; R ) by saying that for z ∈ V (Ω; R ), v k = Θ( z ) is asolution of the problem − Z Ω S ( ∇ h v k ) : ∇ h φ d x = < F ( z ) , φ > + < S , φ > for any φ ∈ V (Ω; R ) , (5.2)where < S , φ > = Z Ω S ( ∇ h u B ) : ∇ h φ d x, < F ( z ) , φ > = Z Ω D t ( η ( z ) b v k ) · b φ d x + X K ∈T X σ ∈E ( K ) ∩E int Z σ F σ,K ( η ( z ) b v k , u k ) · b φ d S x + X K ∈T X σ ∈E ( K ) Z σ η ( z ) b u k · n σ,K u B · b φ d S x + X σ ∈E out Z σ η ( z ) u B,σ · n σ b v k · b φ d S x + X σ ∈E in Z σ ̺ B u B,σ · n σ b v k · b φ d S x − Z Ω p h ( η ( z ))div h φ d x + κh ω X σ ∈E int Z σ [[ η ( z k )]] σ { b v k } σ [[ b φ ]] σ d S x . Clearly, u k = v k + u B , ̺ k = η ( v k ), with v k any fixed point of the map Θ solves the algebraic problem(3.1–3.4).We apply to the map Θ Schaeffer’s fixed point theorem (see Lemma 4.1) with Z = V (Ω; R ). We haveonly to verify that the set { z ∈ Z | z solves (5.2) with r.h.s. Λ F ( z ) + Λ S , Λ ∈ [0 , } is non empty (indeed 0 belongs to the set) and is bounded. The boundedness, follows from the energybalance. Indeed, the same calculation –as we effectuate later in Section 6 to derive the energy estimates forthe numerical scheme– performed on the modified problem (5.2) with v k = z , yields (6.4) with v k replacedby z , where all terms are multiplied by Λ except the term ∆ t R Ω S ( ∇ h z ) : ∇ h z d x at its left hand side. Thisreadily implies the boundedness of the latter set by the argumentation in the proof of Lemma 6.1.This completes the proof of Theorem 3.1.
1. We take in the continuity equation (3.3) the test function φ = 1 and obtain the conservation of mass(see equation (6.3)).2. We shall mimick the derivation of the energy iequality in the contiuous case (cf. [52]) and compute(3.4) φ = v k + (3.3) φ = − | b v k | + (3.3) φ = H ′ h ( ̺ k ) , (6.1)where H h ( ̺ ) = H ( ̺ ) + ˜ κh η ̺ . (6.2)Regrouping conveniently various terms, we get the following identities:(a) Contribution of the terms with the time derivatives to (6.1): Z Ω (cid:16) D t [ ̺ k v k ] · b v k − D t [ ̺ k ] | b v k | + D t [ ̺ k ] H ′ h ( ̺ k ) (cid:17) d x = Z Ω D t h ̺ k | b v k | + H h ( ̺ k ) i d x +∆ t Z Ω ̺ k − | D t [ b v k ] | d x + 1∆ t Z Ω (cid:16) H h ( ̺ k − ) − H ′ h ( ̺ k )( ̺ k − − ̺ k ) − H h ( ̺ k ) (cid:17) d x. φ = v k + (3.3) φ = − | b v k | : X K ∈T X σ ∈E ( K ) ∩E int Z σ F σ,K ( ̺ k b v k , u k ) · b v k d S x − X K ∈T X σ ∈E ( K ) ∩E int Z σ F ( ̺ k , u k ) | b v k | d S x = 12 X σ ∈E int Z σ | Up σ ( ̺ k , u k ) | [[ b v k ]] σ . (c) Contribution of the artificial density dissipation to (3.4) φ = v k + (3.3) φ = − | b v k | : κh ω X σ ∈E int Z σ (cid:16) [[ ̺ k ]] σ { b v k } σ · [[ b v k ]] σ −
12 [[ ̺ k ]] σ [[ | b v k | ]] σ (cid:17) d S x = 0 . (d) Contribution of the pressure term in the momentum and of the convective term in the continuityequation to (3.4) φ = v k + (3.3) φ = H ′ ( ̺ k ) : − Z Ω p h ( ̺ k )div h v k d x + X K ∈T X σ ∈E ( K ) ∩E int Z σ F σ,K ( ̺ k , u k ) H ′ h ( ̺ k )d S x + X σ ∈E out Z σ ̺ k u B,σ · n σ H ′ h ( ̺ k )d S x + X σ ∈E in Z σ ̺ B,σ u B,σ · n σ H ′ h ( ̺ k )d S x = X K ∈T X σ ∈E ( K ) Z σ p h ( ̺ k ) u B,σ · n σ,K d S x + X σ ∈E int Z σ (cid:16) H h ( ̺ k, − ) − H ′ h ( ̺ k, + )( ̺ k, − − ̺ k, + ) − H h ( ̺ k, + ) (cid:17) | u kσ · n σ | d S x + X σ ∈E out Z σ H h ( ̺ k ) u B,σ · n σ d S x + X σ ∈E in Z σ H h ( ̺ B ) u B,σ · n σ d S x + X σ ∈E in Z σ (cid:16) H h ( ̺ B ) − H ′ ( ̺ k )( ̺ B − ̺ k ) − H h ( ̺ k ) (cid:17) | u B,σ · n σ | d S x . (e) Contribution of the artificial viscosity term in the continuity equation to to (3.3) φ = H ′ ( ̺ k ) is κh ω X σ ∈E int Z σ [[ ̺ k ]] σ [[ H ′ h ( ̺ k )]] σ d S x . Putting together items 2(a)–2(e) while taking into account the evident contribution of the dissipative termsfor the velocity in the momentum equation to (3.4) φ = v k , we get the energy balance. We gather the abovecalculations in the following lemma. Lemma 6.1 [Mass conservation and energy balance]
Suppose that the pressure satisfies assumptions (1.6).Then any solution of ( ̺, u ) of the algebraic system (3.1–3.4) satisfies for all m = 1 , . . . , N the following: . Mass conservation ̺ m > , Z Ω ̺ m d x + ∆ t m X k =1 X σ ∈E out Z σ ̺ k u B,σ · n σ d S x = Z Ω ̺ d x − ∆ t m X k =1 X σ ∈E in Z σ ̺ B u B,σ · n σ d S x ; (6.3)
2. Energy balance. There exist ̺ k − ,k ∈ Q (Ω) , ̺ k − ,kK ∈ [min { ̺ k − K , ̺ kK } , max { ̺ k − K , ̺ kK } ] , K ∈ T and ̺ k,σ ∈ [min { ̺ − σ , ̺ + σ } , max { ̺ − σ , ̺ + σ } ] , σ ∈ E int , k = 1 , . . . , N , such that Z Ω ̺ m | b v m | d x + Z Ω H h ( ̺ m ) d x + ∆ t m X k =1 Z Ω S ( ∇ h v k ) : ∇ h v k d x (6.4)+ 12 m X k =1 Z Ω (cid:16) ̺ k − | b v k − b v k − | + H ′′ h ( ̺ k − ,k ) | ̺ k − ̺ k − | (cid:17) d x + ∆ t m X k =1 X σ ∈E int Z σ | Up σ ( ̺ k , u k ) | [[ b v k ]] σ d S x + κh ω ∆ t m X k =1 X σ ∈E int Z σ [[ ̺ k ]] σ [[ H ′ h ( ̺ k )]] σ d S x + ∆ t m X k =1 X σ ∈E int Z σ H ′′ h ( ̺ k,σ )[[ ̺ k ]] σ | u kσ · n σ | d S x + ∆ t m X k =1 X σ ∈E out Z σ ̺ k u kσ · n σ | b v k | d S x +∆ t m X k =1 X σ ∈E out Z σ H h ( ̺ k ) u B,σ · n σ d S x + ∆ t m X k =1 X σ ∈E in E H h ( ̺ B | ̺ k ) | u B,σ · n σ | d S x = Z Ω h ̺ | b v | + H h ( ̺ ) i d x + ∆ t m X k =1 X σ ∈E in Z σ H h ( ̺ B ) | u B,σ · n σ | d S x − ∆ t m X k =1 Z Ω S ( ∇ h u B ) : ∇ h v k d x − ∆ t m X k =1 Z Ω p h ( ̺ k )div h u B d x − ∆ t m X k =1 Z Ω ̺ k b u k · ∇ h u B · b v k d x − ∆ t m X k =1 X σ ∈E in Z σ ̺ B u B,σ · n σ | b v k | d S x . We are now in position to deduce from the energy balance (6.4) the uniform estimates. This will bedone by using the Gronwall lemma. Before starting, it would be convenient to mention some properties ofthe pressure p and of the Helmholtz function H which can be deduced from assumptions (1.6–1.8). Remark 6.1
1. We know that H is strictly convex (since p ′ ( ̺ ) > ). Consequently p is strictly convexas well (since ap − H is convex). We may thus suppose, without loss of generality, that both H − ap and ap − H are strictly convex functions.2. It is easy to deduce from the convexity of function ap − H that ∃ R > , ∀ ̺ > R, ̺ γ − < ∼ H ′ ( ̺ ) , ̺ γ − < ∼ p ′ ( ̺ ) where γ = 1 + 1 /a, (6.5) see [1, formula (2.3)]. This observation explains the value of γ in Theorem 3.2. . Convexity of ap − H and p (0) = H (0) = 0 yield ∀ ̺ ∈ (0 , ∞ ) , p ( ̺ ) < ∼ a H ( ̺ ) + c̺, c = ap ′ (0) − H ′ (0) a . (6.6) This observation allows to prove later Lemma 10.2 and consequently the error estimates.4. Since H − ap is convex, there are numbers d ≥ , e ≥ such that inf ̺> ( H ( ̺ ) − ap ( ̺ )) ≥ , lim ̺ →∞ ( H ( ̺ ) − ap ( ̺ )) = ∞ , where H ( ̺ ) = H ( ̺ ) + d̺ + e. (6.7) This observation allows to derive the uniform estimates from the energy balance (6.4).
With this observation at hand, for the purpose of obtaining the uniform estimates, we can replace in theenergy balance at its left hand side Z Ω H h ( ̺ m ) d x by Z Ω H h ( ̺ m ) d x, H h ( ̺ ) = H ( ̺ ) + ˜ κh η ̺ and add the non negative term ∆ t m X k =1 X σ ∈E out Z σ d̺ k u B,σ · n σ d S x . Likewise we can replace at its right hand side Z Ω H h ( ̺ ) d x by Z Ω H h ( ̺ ) d x and add the term − ∆ t m X k =1 X σ ∈E in Z σ d̺ B u B,σ · n σ d S x . Due to the balance of mass (6.3) the energy balance modified in this way is again satisfied as an identity.Now, we shall write the term R Ω H h ( ̺ m ) d x at the left hand side of the modified energy balance as asum Z Ω ( H h ( ̺ m ) − ap ( ̺ m )) d x + a Z Ω p ( ̺ m ) d x In view of this, in order to deduce from the modified energy balance the uniform estimates, we want to finda bound from above of its right hand side by the expression1 δ Z τ m Z Ω (cid:16) ̺ | b v | + p ( ̺ ) + ˜ κh η ̺ (cid:17) d x d t + δ Z τ m Z Ω |∇ h v | d x d t + C, τ m = m ∆ t with some C > T and Ω and on δ >
0, with the goal to employ the Gronwallinequality.The first line contains only the data. The third term can be treated by the Cauchy-Schwartz and Younginequalities. The fourth term is easy as well. For the fifth term, we have, Z τ m Z Ω ̺ b u · ∇ h u B · b v d x d t < ∼ τ m + Z τ m Z Ω ̺ b v d x d t. | X σ ∈E in Z σ ̺ B u B,σ · n σ | b v | d S x | < ∼ X σ ∈E in k b v k L ( σ ) (6.8) < ∼ h − X K ∈T , K ∩E in = ∅ k v k L ( K ) < ∼ h X K ∈T , K ∩E in = ∅ k∇ x v k L ( K ) < ∼ h k∇ h v k L (Ω) . Indeed the first inequality is the Cauchy-Schwartz inequality for the integrals over σ . The second inequalityuses the discrete trace estimate (4.15) and the Jensen’s inequality (4.1). Finally, we use one of the inequalities(4.3) in combination with the fact that v σ = 0 for any σ ∈ E in . In the following lemma we shall gather all estimates which can be deduced from Lemma 6.1.
Lemma 6.2
We denote I = [0 , T ] , E = sup h ∈ (0 , E ,h , E ,h = Z Ω (cid:20) ̺ h | b u h | + H h ( ̺ h ) (cid:21) d x. We suppose that the pressure satisfies assumption (1.6–1.8) and that h = ∆ t .Then there exists a number d , d := d ( k u B k C (Ω) , k r B k C (Ω) , E , T, Ω) < ∼ , such that any solution ( ̺, u ) =( ̺ h, ∆ t , u h, ∆ t ) of the discrete problem (3.1–3.4) admits the following estimates: ̺ > , k ̺ k L ∞ ( I ; L γ (Ω)) < ∼ , k p ( ̺ ) k L ∞ (0 ,T ; L (Ω)) < ∼ , k H ( ̺ ) k L ∞ (0 ,T ; L (Ω)) < ∼ , (6.9) k ̺ | u B · n | /γ k L γ ( I ; L γ ( ∂ Ω)) < ∼ , k H ( ̺ ) | u B · n |k L ( I ; L ( ∂ Ω)) < ∼ , (6.10)˜ κh η/ k ̺ k L ∞ ( I ; L (Ω)) < ∼ , (6.11)˜ κh η/ k ̺ | u B · n | / k L ( I ; L ( ∂ Ω)) < ∼ , (6.12)sup τ ∈ (0 ,T ) k√ ̺ b u ( τ, · ) k L ∞ ( I ; L (Ω) < ∼ , (6.13) k∇ h v k L ( I × Ω) < ∼ , k v k L ( I ; L (Ω)) < ∼ , (6.14)˜ κh η Z T Z Ω ∆ t | D t ̺ | d x d t = ˜ κh η m X k =1 Z Ω (cid:12)(cid:12)(cid:12) ̺ k − ̺ k − (cid:12)(cid:12)(cid:12) d x < ∼ , (6.15) Z T Z Ω ∆ t̺ ( t − ∆ t ) | D t v ( t ) | d x d t = m X k =1 Z Ω ̺ k − (cid:12)(cid:12)(cid:12)b v k − b v k − (cid:12)(cid:12)(cid:12) d x < ∼ , (6.16) X σ ∈E int Z T Z σ | Up σ ( ̺, u ) | [[ v ]] σ d S x d t < ∼ , (6.17) Here we use the definition of Π V and the property (2.1) of the mesh. The latter estimate follows from the first one by the discrete Sobolev inequality (4.11). κh η X σ ∈E int Z T Z σ [[ ̺ ]] σ | u σ · n | dS x d t < ∼ , (6.18) κ ˜ κh ω + η X σ ∈E int Z T Z σ [[ ̺ ]] σ dS x d t < ∼ . (6.19) Remark 6.2
Some other estimates can be deduced if we consider particular situations enumerated in Re-mark 3.2. They are essential in the proof of the versions of Theorem 3.2 in those particular situations.Suppose that the pressure satisfies in addition to (1.6–1.8) also (3.15). Then we have:1. If in (3.15), π > and γ ≥ we have in addition to (6.18), X σ ∈E int Z T Z σ [[ ̺ ]] σ | u σ · n | dS x d t < ∼ , (6.20) and, in addition to (6.15) Z T Z Ω ∆ t | D t ̺ | d x d t = m X k =1 Z Ω (cid:12)(cid:12)(cid:12) ̺ k − ̺ k − (cid:12)(cid:12)(cid:12) d x < ∼ . (6.21)
2. If < γ ≤ , we have in addition to (6.18) X σ ∈E int Z T Z σ (cid:16) { ̺ } σ + 1 (cid:17) − γ [[ ̺ ]] σ | u σ · n | dS x d t < ∼ , (6.22) and, in addition to (6.15), m X k =1 Z Ω (cid:16) ̺ k − + ̺ k + 1 (cid:17) − γ (cid:12)(cid:12)(cid:12) ̺ k − ̺ k − (cid:12)(cid:12)(cid:12) d x < ∼ . (6.23)
3. If < γ < and condition (3.17) is satisfied, we have also κh ω X σ ∈E int Z T Z σ [[ ̺ γ/ ]] σ dS x d t < ∼ , (6.24) and consequently κh ω Z T k ̺ k γL γ (Ω) d t < ∼ . (6.25) Having collected all the available uniform bounds, our next task is to verify that our numerical method is consistent with the variational formulation of the original problem. This requires some explanation: We deduce from (3.17) that [[( ̺ k ) γ/ ]] σ < ∼ H ′′ ( ̺ k,σ )( ̺ k, − − ̺ k, + ) and thus the fourth lineof (6.4) yields, in particular, (6.24). The latter implies (6.25) by virtue of Sobolev inequality (4.9) for the broken norms. .1 Continuity equation For φ ∈ C ([0 , T ] × Ω), take b φ ( t ) as a test function in the discrete continuity equation (3.3). We see that Z Ω D t ̺ b φ d x = Z Ω D t ̺φ d x. Further, using the formula (4.22), and noticing that Z Ω ( b φ − φ ) ̺ div h u d x = X K ∈T Z K ( b φ − φ ) ̺ div x u = 0(indeed div h u ( t ) is constant on each element K ), we check without difficulty using formula (4.22) that Z Ω ̺ u · ∇ x φ d x = X K ∈T X σ ∈E ( K ) Z σ ( φ − b φ ) ̺ ( u − u σ ) · n σ,K dS x + Z Ω ( b φ − φ ) ̺ div h u d x − X K ∈T X σ ∈E ( K ) ∩E int Z σ F σ,K ( ̺, u ) b φ d S x + X K ∈T X σ ∈E ( K ) ∩E int Z σ ( b φ − φ )[[ ̺ ]] n σ,K [ u σ · n σ,K ] − d S x + X σ ∈E out Z σ ̺ u B,σ · n σ ( φ − b φ )d S x + X σ ∈E in Z σ ̺ u B,σ · n σ ( φ − b φ )d S x . Consequently, equation (3.3) rewrites Z Ω (cid:16) D t ̺φ − ̺ b u · ∇ x φ (cid:17) d x + Z Γ out ̺ u B · n φ d S x + Z Γ in ̺ B u B · n φ d S x = < R Ch ( t ) , φ > (7.1)in (0 , T ], where < R Ch ( t ) , φ > = X K ∈T X σ ∈E ( K ) ∩E int Z σ ( φ − b φ )[[ ̺ ]] n σ,K [ u σ · n σ,K ] − d S x + X K ∈T X σ ∈E ( K ) Z σ ( b φ − φ ) ̺ ( u − u σ ) · n σ,K dS x + Z Ω ̺ div h u ( φ − b φ ) d x − κh ω X σ ∈E int Z σ [[ ̺ ]] σ [[ b φ ]] σ d S x + X σ ∈E in Z σ ( ̺ − ̺ B ) u B,σ · n σ ( b φ − φ )d S x + Z Ω ̺ ( u − b u ) · ∇ x φ d x. In the above, we have also used the fact that the mesh fits to the inflow-outflow boundary and the definition(2.6) of the projection Π V .Alternatively, in view of (2.13), we can integrate the equation (7.1) over time and rewrite it in the form Z Ω ̺ e φ ( τ ) d x − Z Ω ̺ φ (0) d x − Z τ Z Ω (cid:16) ̺ e ∂ t φ + ̺ b u · ∇ x φ (cid:17) d x d t (7.2)+ Z τ Z Γ out ̺ u B · n φ d S x d t + Z τ Z Γ in ̺ B u B · n φ d S x d t = Z τ < R Ch ( t ) , φ > d t with any τ ∈ (0 , T ] and all φ ∈ C ([0 , T ] × Ω), where ̺ e is defined in (2.13).Our next goal is to estimate conveniently all terms in the remainder < R Ch, ∆ t , φ > . To do this we shalluse the tools evoked in Section 4.1 and employ the bounds (6.9–6.19).26. The first term in < R Ch , φ > is bounded from above by < ∼ h h η X σ ∈E int Z σ [[ ̺ ]] σ | u σ · n σ | dS x i / × h h − η X K ∈T X σ ∈E ( K ) ∩E int Z σ ( b φ − φ ) | u σ · n σ | d S x i / , (7.3)where the first term in the product is controlled by (6.18). As far as for the second term in the product,we have for any γ > h − η (cid:12)(cid:12)(cid:12) X K ∈T X σ ∈E ( K ) ∩E int Z σ ( b φ − φ ) | u σ · n σ | d S x (cid:12)(cid:12)(cid:12) < ∼ h − η X K ∈T X σ ∈E ( K ) ∩E int Z σ k b φ − φ k L ∞ ( σ ) k u k L ( σ ) < ∼ h − η X K ∈T k∇ x φ k L ∞ ( K ) k u k L ( K ) < ∼ h − η k∇ x φ k L ∞ (Ω) k u k L (Ω) , (7.4)where k u k L (Ω) < ∼ k u k L (Ω) is bounded in L ((0 , T )). Indeed, to get the second line we have employedthe H¨older and Jensen inequalities on σ , cf. (4.1), third line employs the trace estimates (4.14– 4.15)and then the H¨older inequality.2. By the same token, the absolute value of the second term in < R Ch , φ > X K ∈T X σ ∈E ( K ) Z σ ( b φ − φ ) ̺ ( u − u σ ) · n σ,K dS x allows for any γ > < ∼ h − η k h η/ ̺ k L (Ω) k∇ h u k L (Ω) k∇ x φ k L ∞ (Ω) < ∼ h − η k∇ h u k L (Ω) k∇ x φ k L ∞ (Ω) , (7.5)where we have used (4.3). The third term and the last term in < R Ch , φ > admit the same estimate.3. The artificial density diffusion term in < R Ch , φ > can be processed in the similar way, < ∼ κh ω k∇ x φ k L ∞ (Ω) X σ ∈E int k ̺ k L γ ( σ ) | σ | /γ ′ < ∼ κh ω k ̺ k L γ (Ω)) k∇ x φ k L ∞ (Ω) (7.6)for any value γ > The boundary term in < R Ch ( t ) , φ > is controlled in the same way (without loss of generality, weperform the calculation just for the first of them, and with ̺ only -instead of with ( ̺ − ̺ B )): < ∼ h k∇ φ k L ∞ (Ω) X K ∈T , E ( K ) ∩E in = ∅ X σ ∈E in k ̺ k L ( σ ) < ∼ k∇ φ k L ∞ (Ω) k ̺ k L ( U ) |U | / < ∼ h / − η/ k∇ φ k L ∞ (Ω) , (7.7)where U = ∪ K ∈T , E ( K ) ∩E in = ∅ K . Indeed, in the passage from the first to the second line, we have usedtrace estimates (4.15), and for the rest the H¨older inequalities as well as the fact that the Lebesguemeasure |U | of U is ≈ h .5. At one place of the convergence proof, we shall need a slightly improved version of the above esti-mates with test functions φ ∈ L p (0 , T ; W ,p (Ω)) with p ”large”. To achieve this, we observe thatthe final bound in (7.4) can be replaced by < ∼ h − η k∇ x φ ( t ) k L p (Ω) k u ( t ) k L (Ω) provided p ≥ /
5, fi-nal bound in (7.5) by < ∼ h p − p − η k h η/ ̺ k L (Ω) k u ( t ) k L (Ω) k∇ x φ ( t ) k L p (Ω) provided p − p − η > < ∼ κh ω k ̺ ( t ) k L γ (Ω)) k∇ x φ ( t ) k L p (Ω) provided γ + p > Lemma 7.1 [Consistency for the continuity equation]
Let h = ∆ t , ˜ κ = 1 and η ∈ (0 , . Let the pressuresatisfy the hypotheses (1.6–1.8).1. Then any numerical solution of problem (3.1–3.4) satisfies the continuity equation in the variationalform (7.2) with any φ ∈ C ([0 , T ] × Ω) , where (cid:12)(cid:12)(cid:10) R Ch , φ (cid:11)(cid:12)(cid:12) < ∼ h α C r Ch k φ k L ∞ (0 ,T ; W , ∞ (Ω)) , k r Ch k L / ((0 ,T )) < ∼ with α C = − η > if κ = 0 and α C = min { − η , ω } > if κ = 1 .2. There is p ( η ) > such that for all p > p ( η ) Z τ (cid:12)(cid:12)(cid:10) R Ch , φ (cid:11)(cid:12)(cid:12) d t < ∼ h α C k φ k L p (0 ,T ; W ,p (Ω)) , (7.9) with some α C = α C ( p ) > and all φ ∈ C c ((0 , T ) × Ω) . Remark 7.1
Let ˜ κ = 0 and let p satisfy in addition to (1.6–1.8) also condition (3.15).1. Formulas (7.8–7.9) with α C > and α C > remain still valid under condition / < γ < providedalso (3.17) is satisfied. In this case α C = min n , γ − γ , γ ′ o if κ = 0 and it is minimum of the above value and ω if κ = 1 . Also, there is p = p ( γ ) > such that (7.9)holds for all p > p with some α C = α C ( p ) > .Indeed, we can use estimates (6.22) instead of (6.18) when calculating (7.3) which now yields, < ∼ h X σ ∈E int Z σ (cid:16) { ̺ } σ + 1 (cid:17) − γ [[ ̺ ]] σ | u σ · n σ | dS x i / × h X K ∈T X σ ∈E ( K ) ∩E int Z σ (cid:16) { ̺ } σ + 1 (cid:17) − γ ( b φ − φ ) | u σ · n σ | d S x i / and the estimate (we use systematically the ”negative” interpolation (4.17)), < ∼ h k ̺ k − γ L γ (Ω) k∇ φ k L ∞ (Ω) k u k / L (Ω) A h < ∼ h B h k∇ φ k L ∞ (Ω) , B h bounded in L / ((0 , T )) for the first term in R Ch .Further, we use estimate k ̺ k L (Ω) < ∼ h min { , − γ } k ̺ k L γ (Ω) , when evaluating term corresponding to (7.5), which gives the bound < ∼ h { , − γ } k∇ φ k L ∞ (Ω) A h , A h bounded in L ((0 , T ))28 n these cases. Estimate (7.6) of the artificial diffusion term remains in force, while the the estimate of the boundaryterm (7.7) can be replaced by < ∼ k∇ φ k L ∞ (Ω) k ̺ k L γ ( U ) |U | /γ ′ < ∼ h /γ ′ k∇ φ k L ∞ (Ω) .
2. Similar reasoning can be carried out if γ ≥ and if in (3.15), π > . In this case, in (7.8), α C = min n , γ ′ o if κ = 0 and it is minimum of the above value and ω if κ = 1 . The next step is to take ˜ φ ( t ) , φ ∈ C c ([0 , T ] × Ω; R ) , as a test function in the discrete momentum equation (3.4). Seeing that, in accordance with (4.12), (4.13), Z Ω S ( ∇ h v ) : ∇ h ˜ φ d x = Z Ω S ( ∇ h v ) : ∇ x φ d x, Z Ω p h ( ̺ )div h ˜ φ d x = Z Ω p h ( ̺ )div x φ d x, we may rewrite (3.4)–by using the formula (4.22) and rearanging conveniently several terms– in the form Z Ω D t ( ̺ b v ) · φ d x − Z Ω ̺ b u ⊗ b v : ∇ x φ d x − Z Ω ̺ b u · ∇ x φ · b u B d x + Z Ω ̺ b u · ∇ h ( u B · φ ) d x (7.10)+ Z Ω S ( ∇ h u ) : ∇ x φ d x − Z Ω p h ( ̺ )div x φ d x = < R Mh, ∆ t , φ > in (0 , T ] , where < R Mh, ∆ t , φ > = Z Ω D t ( ̺ b v ) · ( φ − b ˜ φ ) d x + X K ∈T X σ ∈E ( K ) ∩E int Z σ ( φ − b ˜ φ ) · [[ ̺ b v ]] σ, n σ,K [ u σ · n σ,K ] − d S x + X K ∈T X σ ∈E ( K ) Z σ ̺ ( b ˜ φ − φ ) · b v ( u − u σ ) · n σ,K d S x + Z Ω ̺ ( φ − b ˜ φ ) · b v div h u d x + Z Ω ̺ ( u − b u ) · ∇ x φ · b v d x + Z Ω ̺ b u · ∇ h u B · ( φ − b ˜ φ ) d x + Z Ω ̺ b u · ∇ φ · ( u B − b u B ) d x − κh ω X σ ∈E int Z σ [[ ̺ k ]] σ { b u k } σ [[ b ˜ φ ]] σ d S x + X σ ∈E in Z σ ( ̺ B − ̺ ) u B,σ · n σ b v · ( φ − b ˜ φ ) = X i =1 I i Alternatively, in view of (2.13), equation (7.10) can be rewritten as follows, Z Ω ̺ b v f · φ ( τ, x ) d x − Z Ω ̺ b v · φ (0 , x ) d x − Z τ Z Ω h ̺ b v f · ∂ t φ + (cid:16) ̺ b u ⊗ b u + p h ( ̺ ) I (cid:17) : ∇ x φ i d x d t (7.11) This is the place when γ > / h positive. Here the situation could be still handledfor 1 < γ ≤ / κ = 1 by using estimate (6.25). Similar terms in the remainder of the consistency formulation of themomentum equation, however, impose the restriction γ > / κ = 1. It is therefore useless to treat this situation here. Z τ Z Ω ̺ b u · ∇ h ( u B · φ ) d x + Z τ Z Ω S ( ∇ h u ) : ∇ x φ d x d t = Z τ < R Mh, ∆ t ( t ) , φ > d t with any τ ∈ (0 , T ] and all φ ∈ C c ([0 , T ] × Ω , R ).Our goal is to estimate conveniently all terms in < R Mh, ∆ t , φ > . As in Section 7.1, we shall use the toolsreported in Section 4.1 and uniform estimates for numerical solutions derived in Lemma 6.2.1. Most terms in < R
Mh, ∆ t , φ > contain the expression b ˜ φ − φ = Π Q ( ˜ φ − φ ) + b φ − φ . We notice that byvirtue of (4.3), (4.5), (4.6), k φ − b ˜ φ k L q (Ω) < ∼ h k∇ x φ k L q (Ω) , ≤ q ≤ ∞ . Estimate of the first term in R Mh,h (term I ) : For the error in the time derivative, we obtain (cid:12)(cid:12)(cid:12)(cid:12)Z Ω D t ( ̺ b v ) · ( φ − b ˜ φ ) d x (cid:12)(cid:12)(cid:12)(cid:12) < ∼ √ h (∆ t ) − / A h , A h bounded in L (0 , T ) . Indeed, we split the term R Ω D t ( ̺ k b v k ) · ( φ − b ˜ φ ) d x into two parts, Z Ω p ̺ k − p ̺ k − v k − v k − ∆ t · ( φ − b ˜ φ ) d x + Z Ω ̺ k − ̺ k − ∆ t v k · ( φ − b ˜ φ ) d x, (7.12)where, by virtue of H¨older’s inequality, for any γ >
1, the first term is bounded by < ∼ h (∆ t ) − / k ̺ k − k / L γ (Ω) ∆ t Z Ω ̺ k (cid:18) v k − − v k − ∆ t (cid:19) d x ! / k∇ x φ k L γγ − (Ω) (7.13)while the second one is < ∼ h − η/ (∆ t ) − / ∆ t h η Z Ω (cid:18) ̺ k − ̺ k − ∆ t (cid:19) d x ! / k v k k L (Ω; R ) k∇ x φ k L (Ω) , (7.14)where the first integrals of both expressions are controlled by means of (6.15) and (6.16).3. Estimate of the second term in R Mh,h (term I ) : This essentially amounts to estimate two terms, I , = X K ∈T X σ ∈E ( K ) ∩E int Z σ ̺ + σ, n σ,K ( φ − b ˜ φ ) · [[ b v ]] σ, n σ,K [ u σ · n σ,K ] − d S x and I , = X K ∈T X σ ∈E ( K ) ∩E int Z σ ( φ − b ˜ φ ) · b v − σ, n σ,K [[ b ̺ ]] σ, n σ,K [ u σ · n σ,K ] − d S x , where | I , | < ∼ h X K ∈T X σ ∈E ( K ) ∩E int Z σ ̺ + σ, n σ,K [[ b v ]] σ, n σ,K | [ u σ · n σ,K ] − | d S x i / × h X K ∈T X σ ∈E ( K ) ∩E int Z σ ̺ + σ, n σ,K ( φ − b ˜ φ ) | [ u σ · n σ,K ] − | d S x i / (7.15)30ith the first term in the product controlled by (6.17), i.e. belonging to L ((0 , T )), while | I , | < ∼ h h η X K ∈T X σ ∈E ( K ) ∩E int Z σ [[ b ̺ ]] σ, n σ,K | [ u σ · n σ,K ] − | d S x i / × h h − η X K ∈T X σ ∈E ( K ) ∩E int Z σ ( φ − b ˜ φ ) · | b v − σ, n σ,K | | [ u σ · n σ,K ] − | d S x i / (7.16)with the first term in the product controlled by (6.18), i.e. belonging also to L ((0 , T )). It will betherefore enough to estimate the expressions under the second square roots of I , , I , , respectively.In fact, it is enough to estimate only the ”leading terms” of these expressions, where we replace u by v .We have h − η (cid:12)(cid:12)(cid:12) X K ∈T X σ ∈E ( K ) ∩E int Z σ ̺ + σ, n σ,K ( φ − b ˜ φ ) | [ v σ · n σ,K ] − | d S x (cid:12)(cid:12)(cid:12) < ∼ h − η k∇ φ k L ∞ (Ω)) X K ∈T X σ ∈E ( K ) ∩E int k ̺ k L ( σ ) k v k L ( σ ) (7.17) < ∼ h − η k∇ φ k L ∞ (Ω)) X K ∈T k ̺ k L ( K ) k v k L ( K ) < ∼ h − η k∇ φ k L ∞ (Ω)) k h η/ ̺ k L (Ω) k v k L (Ω) . Likewise (cid:12)(cid:12)(cid:12) X K ∈T X σ ∈E ( K ) ∩E int Z σ ( φ − b ˜ φ ) | b v − σ, n σ,K | | [ v σ · n σ,K ] − | d S x (cid:12)(cid:12)(cid:12) < ∼ h k∇ x φ k L ∞ (Ω) k v k L (Ω) < ∼ h (∆ t ) − / k (∆ t ) / v k L (Ω) k∇ x φ k L ∞ (Ω) , (7.18)where k (∆ t ) / v k L (Ω) is bounded in L ((0 , T )) by virtue of (4.19).4. The upper bound of the third term in R Mh,h (term I ) is determined by the upper bound of | X K ∈T X σ ∈E ( K ) Z σ ̺ ( φ − b ˜ φ ) · b v ( v − v σ ) · n σ,K d S x | , which is < ∼ h k ̺ k L (Ω) k v k L ∞ (Ω) k∇ h v k L (Ω) k∇ x φ k L ∞ (Ω) (7.19) < ∼ h − η k h η/ ̺ k L (Ω) k v k L (Ω) k∇ h v k L (Ω) k∇ x φ k L ∞ (Ω) , where k v k L (Ω) k∇ h v k L (Ω) is bounded in L ((0 , T )).5. The bounds of term I in R Mh,h are determined by the bounds of (cid:12)(cid:12)(cid:12) Z Ω ̺ ( φ − b ˜ φ ) · b v div h v d x (cid:12)(cid:12)(cid:12) . They are exactly the same as in the previous case. The same is true for the terms I – I , since theyhave the same structure. 31. The the artificial viscosity term | I | = κh ω | P σ ∈E int R σ [[ ̺ k ]] σ { b v k } σ [[ b φ ]] σ d S x | is bounded by < ∼ κh ω k ̺ k L (Ω) k v k L (Ω) k∇ φ k L ∞ (Ω) < ∼ κh ω − η/ k h η/ ̺ k L (Ω) k v k L (Ω) k∇ φ k L ∞ (Ω) (7.20)where k v k L (Ω) < ∼ k v k L (Ω) is bounded in L ((0 , T )).7. The last term in R Mh,h to be evaluated is the boundary term ( I ) whose decay with h is determined by X σ ∈E in Z σ ̺ u B,σ · n σ b v · ( b ˜ φ − φ )d S x . We have, by the same reasoning as in (7.7), the bound < ∼ h k ̺ k L (Ω) k∇ h v k L (Ω) kk∇ x φ k L ∞ (Ω) < ∼ h − η/ k h η/ ̺ k L (Ω) |∇ h v k L (Ω) k∇ x φ k L ∞ (Ω) . (7.21) Lemma 7.2 [Consistency for the momentum equation]
Let h = ∆ t , ˜ κ = 1 and η ∈ (0 , / . Let the pressuresatisfy the hypotheses (1.6–1.8). Then any numerical solution of problem (3.1–3.4) satisfies the momentumequation in the variational form (7.11) with any φ ∈ C ([0 , T ] × Ω) , where (cid:12)(cid:12)(cid:10) R Mh,h ( t ) , φ (cid:11)(cid:12)(cid:12) < ∼ h α M r Mh k φ k L ∞ (0 ,T ; W , ∞ (Ω)) , k r Mh k L ((0 ,T )) < ∼ with α M = − η > if κ = 0 and α M = min { ω − η , − η } if κ = 0 , η < ω . Remark 7.2
Let ˜ κ = 0 and let p satisfy in addition to (1.6–1.8) also condition (3.15).1. Formula (7.22) with α M > remains still valid under condition / < γ < provided (3.17) issatisfied. In this case, if κ = 0 and γ > / , α M = min n , γ − γ , γ − γ o and if κ = 1 , / < γ ≤ / and ω ∈ (0 , γ − , α M = min n , γ − γ , − γ − ω γ , ω o . Here are some hints (we suppose ∆ t = h and we use systematically the ”negative interpolation” (4.17),(4.19)):(a) When estimating the second term in formula (7.12), the control (6.15) must be replaced by thecontrol (6.23). Under condition h = ∆ t , this amounts to replace (7.14) by < ∼ h − " ∆ t Z Ω (cid:16) ̺ k + ̺ k − + 1 (cid:17) − γ (cid:18) ̺ k − ̺ k − ∆ t (cid:19) d x / h Z Ω ( ̺ k + ̺ k − +1) − γ | v k | | b ˜ φ − φ | d x i / < ∼ h k ̺ k − γγ L γ (Ω) A h k v k k L (Ω) k∇ x φ k L ∞ (Ω) , provided γ ≥ / , where A h k v k k L (Ω) is bounded in L ((0 , T )) . b) Decomposition (7.15) remains unchanged and yields the bound < ∼ h k ̺ k L γ (Ω) k u k L (Ω) A h k∇ x φ k L ∞ (Ω) where A h k u k L (Ω) is bounded in L ((0 , T )) , cf. (6.14), (6.17), and (7.18).Decomposition (7.16) takes the form, h X K ∈T X σ ∈E ( K ) ∩E int Z σ (cid:16) { ̺ } σ + 1 (cid:17) − γ [[ b ̺ ]] σ, n σ,K | [ u σ · n σ,K ] − | d S x i / × h X K ∈T X σ ∈E ( K ) ∩E int Z σ ( { ̺ } σ + 1) − γ ( φ − b ˜ φ ) · | b v − σ, n σ,K | | [ u σ · n σ,K ] − | d S x i / , where the first term is controlled by (6.22). Under assumption ∆ t = h , the expression under thesecond square root (where we take v on place of u without loss of generality) admits the bound < ∼ h +min { , pm − } k ̺ k − γL γ (Ω) k (∆ t ) / v k L (Ω) k∇ x φ k L ∞ (Ω) , where p m + − γγ = 1 and k (∆ t ) / v k L (Ω) is bounded in L ((0 , T )) and the power of h is positiveprovided γ > / . (c) The estimate (7.19) is modified as follows, < ∼ h k ̺ k / L ∞ (Ω) k ̺ b v k / L (Ω) k∇ h v k L (Ω) k∇ x φ k L ∞ (Ω) < ∼ h − γ k∇ h v k L (Ω) k∇ x φ k L ∞ (Ω) , where k∇ h v k L (Ω) is bounded in L ((0 , T )) and the power of h is positive provided γ > / .Alternatively, if κ = 1 , we have in this case the bound < ∼ h − γ − ω γ (cid:16) h ωγ k ̺ k L γ (Ω) (cid:17) / k ̺ b v k / L (Ω) k∇ h v k L (Ω) k∇ x φ k L ∞ (Ω) , where (cid:16) h ωγ k ̺ k L γ (Ω) (cid:17) / k ̺ b v k / L (Ω) k∇ h v k L (Ω) is bounded in L ((0 , T )) and the power of h ispositive provided ω ∈ (0 , γ − . In the above, we have used, in particular, estimate (6.25).(d) The estimate (7.20) now reads < ∼ κh ω k ̺ k L γ (Ω) k v k L (Ω) k∇ φ k L ∞ (Ω) , provided γ ≥ / .(e) Finally, by the argumentation of (7.7), estimate (7.21) can be modified as < ∼ h γ − γ k ̺ k L γ (Ω) k v k L (Ω) k∇ x φ k L ∞ (Ω) .
2. Similar reasoning can be carried out if γ ≥ and in (3.15), π > . In this case, we can take alsoadvantage of estimates (6.20) and (6.21). We obtain, at least, α M = 12 if κ = 0 and α M = min { , ω } if κ = 1 . This is the place, where we cannot get a positive power of h for γ ≤ / κ = 1 and employ estimate (6.25).This is not due to the boundary conditions. Indeed, the same limitation is observed for the no-slip boundary conditions. .3 Energy balance Neglecting several positive terms at the left hand side and using the definition of Π V and (2.1) in order toreplace u B,σ = ˜ u B,σ by u B in the remaining boundary integrals, we can rewrite the energy identity (6.4) inthe following form, h Z Ω (cid:16) ̺ | b v | + H h ( ̺ ) (cid:17) d x i τ + Z τ Z Ω S ( ∇ h u ) : ∇ h v d x d t (7.23)+ Z τ Z Γ out H h ( ̺ ) u B · n d S x d t < ∼ − Z τ Z Γ in H h ( ̺ B ) u B · n d S x d t − Z τ Z Ω (cid:16) ̺ b u ⊗ b u + p h ( ̺ ) I (cid:17) : ∇ h u B d x d t + Z τ Z Ω ̺ b u · ∇ h u B · b u B d x d t + R Eh [ ̺, u ] , where h R Eh [ ̺, u ] i ( τ ) = Z τ m τ h Z Γ in H h ( ̺ B ) | u B · n | d S x + Z Γ out ( d̺ + e ) | u B · n | d S x − Z Ω (cid:16) ̺ b u ⊗ b u + p h ( ̺ ) I (cid:17) : ∇ h u B d x + Z Ω ̺ b u · ∇ h u B · b u B d x i d t, τ ∈ ( τ m − , τ m ] , m = 1 , . . . , N ;whence (cid:12)(cid:12)(cid:12) R Eh [ ̺, u ] (cid:12)(cid:12)(cid:12) < ∼ h, where ∆ t = h (7.24)by virtue of (6.9–6.14), cf. also (6.3), (6.8), (6.7). We denote by [ ̺ h , u h ], h > r , u ], and there is a positively semi-definite tensor measure R (which we shall construct as well), such that the couple [ r , u ] and the associated R is a dissipative solutionof the continuous problem (1.1–1.5) in the sense of Definition 1.1. Recalling regularity (1.3–1.4) of the initial and boundary data, we deduce from (4.2–4.3), (4.5) (4.6), ̺ B,h → r B in L q (Ω) and L q ( ∂ Ω), 1 ≤ q ≤ ∞ , (8.1) b u B,h → u B , u B,h → u B , ∇ h u B,h → ∇ x u B in L q (Ω), 1 ≤ q ≤ ∞ ,̺ h → r in L q (Ω) , u h → u in L q (Ω) , ≤ q ≤ ∞ . By virtue of estimate (6.14) and Jensen’s inequlaity (4.1), k b v h k L (0 ,T ; L (Ω)) < ∼ . (8.2)Consequently, b v h ⇀ v , v h ⇀ v in L (0 , T ; L (Ω)) (8.3) All convergences in this section hold for a chosen subsequence of the original sequence; for the sake of simplicity, we do notrelabel. Z Ω v h ∇ x φ d x = − Z Ω ∇ h v h φ d x + I h + J h , φ ∈ C ((0 , T ) × Ω)) , where I h := X K ∈T X σ ∈E ( K ) Z σ v h n σ,K φ d S x = X K ∈T X σ ∈E ( K ) ∩E int Z σ ( v h − v h,σ ) n σ,K ( φ − φ σ )d S x , J h = Z ∂ Ω v h n φ d S x admit the bounds | I h | < ∼ h k∇ h v h k L (0 ,T ; L (Ω)) k φ k L (0 ,T ; W , (Ω)) , | J h | < ∼ h / k∇ h v h k L (0 ,T ; L (Ω)) k φ k L (0 ,T ; W , (Ω)) (8.4)by virtue of the H¨older inequality, trace estimates (4.14–4.15), the first inequality in (4.3), and in the secondestimate also the standard Sobolev imbedding W , (Ω) ֒ → L (Ω), the fact that v h ∈ V (Ω; R ) and that | ∪ K ∩E ext = ∅ K | < ∼ h. We deduce from this and from (6.14), on one hand, ∇ h v h ⇀ ∇ x v in L (0 , T ; L (Ω)) , and v ∈ L (0 , T ; W , (Ω)), (8.5)and on the other hand, in particular, k∇ x v h k L (0 ,T ; W − , (Ω)) < ∼ , k∇ x b v h k L (0 ,T ; W − , (Ω)) < ∼ , (8.6)where the very latter bound is derived from the previous one by virtue of (4.5).According to (6.9), (6.13), k ̺ h b v h k L ∞ (0 ,T ; L γγ +1 (Ω)) < ∼ . (8.7)It is the consequence of (2.13) and (6.9) that ̺ h f is bounded in L ∞ (0 , T ; L γ (Ω)) . (8.8)Coming back with (8.7) and (6.9) to the continuity equation (7.2), we find, by virtue of Lemma 7.1, thatfor all φ ∈ C c (Ω), Z Ω ̺ h f ( t ) φ d x = A φh ( t ) + B φh ( t ) , where t A φh ( t ) is equi-bounded and equi-continuous in C ([0 , T ]) while B φh → C ([0 , T ]). Consequently,by density of C c (Ω) in L γ ′ (Ω) and by the Arzela-Ascoli type argument, we get ̺ h f → r in C weak ([0 , T ]; L γ (Ω)) . (8.9)We have also, by virtue of (6.9), ̺ h ⇀ ∗ r in L ∞ (0 , T ; L γ (Ω)) . (8.10)The limits (8.9), (8.10) are the same, since k ̺ h f − ̺ h k L ((0 ,T ) × Ω) < ∼ (∆ t ) / h − η/ due to (2.13) and (6.15).Due to (8.7), ̺ h b v h ⇀ ∗ s in L ∞ (0 , T ; L γγ +1 (Ω)) , s = rv (8.11)35here the latter identity follows from Lemma 4.2, where we set r n ≈ ̺ h , v n ≈ b v h , g n ≈ ̺ h b u h , h n ≈ R Ch .(Indeed, one easy checks by using (7.9), (8.7), (8.3), (8.6), (8.9), that assumptions of the Lemma 4.2 aresatisfied. )Finally, due to estimate (6.14) ̺ h ⇀ r in L γ (0 , T ; L γ ( ∂ Ω , | u B · n | d S x )) . (8.12)At this stage we can pass to the limit in the consistency formulation (7.2) of the continuity equation.Due to Lemmma 7.1 Z τ < R Ch , φ > d t → According to (6.9), (6.11), k p h ( ̺ h ) k L ∞ (0 ,T ; L (Ω)) < ∼ . (8.13)In view of (2.13) and (8.7), ̺ h b v h ] is bounded in L ∞ (0 , T ; L γγ +1 (Ω)). (8.14)Coming back with (6.13–6.14), (8.13) to the momentum equation (7.11), we obtain by the same argu-ments as in (8.9), employing now Lemma 7.2 (instead of Lemma 7.1), q h f → q in C weak ([0 , T ]; L γγ +1 (Ω)) , q h := ̺ h b v h , (8.15)where q = rv a.e. in (0 , T ) × Ω . Indeed, the latter identity can be deduced from (8.11) and from the inequality k ̺ h b v h ] − ̺ h b v h k L ((0 ,T ) × Ω) < ∼ (∆ t ) / h − η/ which follows from (2.13) and (6.15–6.16).Due to (8.1) and (8.15), also, m h g → m in C weak ([0 , T ]; L γγ +1 (Ω)) , m h := ̺ h b u h (8.16)where m = ru a.e. in (0 , T ) × Ω and u = v + u B . (8.17)We introduce E ( r, z ) = z ⊗ z r , if z ∈ R , r >
00 if z ∈ R , r = 0+ ∞ if z ∈ R , r < ,E ( ξ ) ( r, z ) := ξ T E ( r, z ) ξ = | z · ξ | r , if z ∈ R , r >
00 if z ∈ R , r = 0+ ∞ if z ∈ R , r < . This is the only point, where we need ”additional” L p -estimate (7.9).
36n the above ξ ∈ R . We observe that ( r, z ) E ( ξ ) ( r, z ) is a lower semi-continuous convex function on R and that E ( ̺ h , b m h ) = ̺ h b u h ⊗ b u h . whence, recalling estimate (6.13), we obtain, by a sequential version of the Banach-Alaoglu-Bourbaki theo-rem E ( ̺ h , b m h ) ⇀ ∗ E ( r , m ) in L ∞ (0 , T ; M (Ω; R × )) (8.18)while ξ T E ( r , b m ) ξ − E ( ξ ) ( r , m ) ≥ , and E ( ξ ) ( r , m ) ∈ L ∞ (0 , T ; L (Ω)) . (8.19)by virtue of Lemma 4.3. Consequently, seeing (8.17), E ( r , m ) = ru ⊗ u ∈ L ∞ (0 , T ; L (Ω)) . We reserve the same treatment to the sequence p h ( ̺ h ) I . Due to (8.13), and since, by (1.8), (3.5), p h isconvex, we have, p h ( ̺ h ) I ⇀ ∗ p ( r ) I in L ∞ (0 , T ; M (Ω; R × )) , (8.20)where p ( r ) − p ( r ) ≥ , p ( r ) ∈ L ∞ (0 , T ; L (Ω)) . Last but not least, we introduce R := (cid:16) E ( r , m ) + p ( r ) I (cid:17) − (cid:16) E ( r , m ) + p ( r ) I (cid:17) ; (8.21)in view of (8.19–8.20),for all ξ ∈ R , ξ T R ξ ∈ L ∞ (0 , T, M + (Ω)) i.e., R ∈ L ∞ (0 , T ; M + (Ω; R × )) . Finally, due to Lemma 7.2, Z τ < R Mh,h , φ > d t → . Now, we are ready to pass to the limit in the momentum equation (7.11), in order to obtain Z Ω q · φ ( τ, x ) d x − Z Ω r v · φ (0 , x ) d x (8.22) − Z τ Z Ω (cid:16) rv · ∂ t φ + (cid:16) ru ⊗ u + p ( r ) I (cid:17) : ∇ x φ (cid:17) d x d t + Z τ Z Ω ru · ∇ ( u B · φ ) d x d t + Z τ Z Ω S ( ∇ x u ) : ∇ x φ d x d t = Z τ Z Ω ∇ x φ : d R ( t )d t. According to (1.11), Z τ Z Ω ru · ∇ x ( u B · φ ) d x d t = − Z τ Z Ω ru B · ∂ t φ d x d t + h Z Ω ru B · φ d x i τ we obtain from (8.22) the formulation (1.12) of the momentum equation. Here and hereafter f ( r , v , ∇ x v , . . . ) denotes a star-weak limit of the sequence f ( ̺ h , v h , ∇ h v h , . . . ). .3 Limit in the energy balance Due to estimate (6.13)Tr[ E ( ̺ h , ̺ h b v h )] = | ̺ h b v h | ̺ h ⇀ ∗ Tr[ E ( r , q )] := (cid:16) q r (cid:17) in L ∞ (0 , T ; M (Ω)) , where, for almost all τ ∈ (0 , T ), (cid:16) q r (cid:17) ( τ ) ≥ Tr[ E ( r ( τ ) , q ( τ ))] = rv ( τ )by virtue of Lemma 4.3. Likewise, by the same token, due to estimate (6.9) H h ( ̺ h ) ⇀ ∗ H ( r ) in L ∞ (0 , T ; M (Ω)) , H ( r ) ≥ H ( r ) . We can thus introduce E = h (cid:16) q r (cid:17) + H ( r ) i − h rv + H ( r ) i ∈ L ∞ (0 , T ; M + (Ω)) . Clearly, in definition of E , one can replace q by m and v by u . Consequently, due to the structural property(1.8) of the pressure p , 0 ≤ D Tr( R ) ≤ E where D = min { / , a/ } . (8.23)Finally, we use (8.12) in conjonction with Lemma 4.3 to show Z τ Z Γ out H ( r ) u B · n d S x d t ≤ lim inf h → Z τ Z Γ out H ( ̺ h ) u B · n d S x d t. We can now pass to the limit in the consistency formulation (7.23) of the energy balance, in order to getenergy inequality in the form (1.13).This finishes the proof of the first part of Theorem 3.2.
Let [ ̺ h , u h ] be any subsequence (not relabeled) of the sequence [ ̺ h , u h ]. Than it contains a subsequenceconverging weakly (in the sense (3.7–3.10)) to a dissipative solution [ r , u ] with defect R .Let [ r, U ] be a strong solution of the problem (1.1–1.5) emanating from the same initial and boundarydata as the dissipative solution [ r , u ]. Then according to Lemmma 1.1, r = r , u = U and the Reynolds defect R = 0.We have one one hand: Any subsequence of the whole sequence [ ̺ h , u h ] admits a weakly convergentsubsequence – in the sense (3.7–3.9) –converging to [ r, U ]; whence the whole sequence weakly converges to[ r, U ].On the other hand R = 0 which implies immediately the convergence (3.11–3.12) by virtue of (8.21),(8.23).This finishes the proof of the first item in the second part of Theorem 3.2.38 The relative energy is a basic tool for showing explicit convergence rate of the strong convergence of numericalsolutions to a strong solution (provided that the latter exist). Our goal is to evaluate the time evolution of E h (cid:16) ̺, b v (cid:12)(cid:12)(cid:12) r, V (cid:17) = E (cid:16) ̺, b v (cid:12)(cid:12)(cid:12) r, V (cid:17) + h η Z Ω ( ̺ − r ) d x where the relative energy functional E ( · , ·|· , · ) is introduced in (1.9), [ ̺, u ] is a numerical solution (of system(3.1–3.4)) and [ r, V ] are test functions in the class (1.14). To achieve this goal we will be mimicking the proofsfrom the ”continuous case”, see [18], [52], [1] using the consistency formulations (7.2), (7.11) and (7.23) ofthe continuity, momentum and energy balance laws. Before starting, to this end, we have, however, makecompatible the consistency formulations of the continuity and momentum equations with the formulation(7.23) of the energy inequality.In view (2.12–2.13), we deduce from (7.2), Z Ω ̺φ ( τ ) d x − Z Ω ̺ φ (0) d x − Z τ Z Ω (cid:16) ̺∂ t φ + ̺ b u · ∇ x φ (cid:17) d x d t + Z τ Z Γ out ̺ u B · n φ d S x d t (10.1)+ Z τ Z Γ in ̺ B u B · n φ d S x d t = < Π Ch, ∆ t ( τ ) , φ >, τ ∈ (0 , T ] , with all φ ∈ C ([0 , T ] × Ω), where < Π Ch, ∆ t ( τ ) , φ > = Z Ω ̺ ( τ )( φ ( τ ) − φ ( τ m )) d x + Z τ m < R Ch, ∆ t , φ > d t + Z τ m Z Ω ( ̺ e − ̺ ) ∂ t φ d x d t + Z τ m τ h Z Ω (cid:16) ̺∂ t φ + ̺ b u · ∇ x φ (cid:17) d x − Z Γ out ̺ u B · n φ d S x − Z Γ in ̺ B u B · n φ d S x i d t, τ ∈ ( τ m − , τ m ] , m = 1 , . . . , N. We have Z τ m Z Ω | ̺ e − ̺ | | ∂ t φ | d x d t < ∼ ∆ t m X k =1 k ̺ k − ̺ k − k L (Ω) k ∂ t φ k L ∞ ((0 ,T ) × Ω) < ∼ (∆ t ) / (cid:16) m X k =1 k ̺ k − ̺ k − k L (Ω) (cid:17) / < ∼ (∆ t ) / h − η/ by virtue of (2.13), where the last inequality follows from (6.15). Employing in estimating of the other terms(6.11), (6.12), (6.13) and (7.8) we get (cid:12)(cid:12)(cid:12) < Π Ch,h ( τ ) , φ > (cid:12)(cid:12)(cid:12) < ∼ h α C k φ, ∇ x φ, ∂ t φ k L ∞ ((0 ,T ) × Ω) . (10.2)Concerning the momentum equation, we deduce from (7.11) by the same token, Z Ω ̺ b v · φ ( τ, x ) d x − Z Ω ̺ b v · φ (0 , x ) d x (10.3) − Z τ Z Ω h ̺ b v · ∂ t φ + (cid:16) ̺ b u ⊗ b u + p h ( ̺ ) I (cid:17) : ∇ x φ i d x d t + Z τ Z Ω ̺ b u · ∇ h ( u B · φ ) d x d t Z τ Z Ω S ( ∇ h u ) : ∇ x φ d x d t = < Π Mh, ∆ t ( τ ) , φ >, τ ∈ (0 , T ] , with all φ ∈ C c ([0 , T ] × Ω , R ), where < Π Mh, ∆ t ( τ ) , φ > = Z Ω ̺ b v ( τ ) · ( φ ( τ ) − φ ( τ m )) d x + Z τ m < R Mh, ∆ t , φ > d t (10.4)+ Z τ m Z Ω ( ̺ b v f − ̺ b v ) ∂ t φ d x d t + Z τ m τ h Z Ω (cid:16) ( ̺ b u ⊗ b u + p h ( ̺ ) I ) : ∇ x φ − ̺ b u · ∇ h ( u B · φ ) (cid:17) d x + Z Ω S ( ∇ h u ) : ∇ x φ d x i d t ; τ ∈ ( τ m − , τ m ] , τ m = m ∆ t, m = 1 , . . . , N ;whence, using (6.11), (6.12), (6.13), (6.15–6.16) and (7.22), we get (cid:12)(cid:12)(cid:12) < Π Mh,h ( τ ) , φ > (cid:12)(cid:12)(cid:12) < ∼ h α M k φ, ∇ x φ, ∂ t φ k L ∞ ((0 ,T ) × Ω) . (10.5) We calculate, (7.23) + (10.2) φ = | V | − H ′ h ( r ) + (10.5) φ = − V
1. We denote U = V + u B and recall that u = v + ˜ u B .2. The above expression gives the following inequality (cid:20)Z Ω E h (cid:16) ̺, b v (cid:12)(cid:12)(cid:12) r, V (cid:17) d x (cid:21) t = τt =0 + Z τ Z Ω S ( ∇ h u ) : ∇ h ( u − U ) d x d t (10.6)+ Z τ Z Γ out (cid:2) H h ( ̺ ) − H ′ h ( r ) ̺ (cid:3) u B · n d S x d t + Z τ Z Γ in (cid:2) H h ( ̺ B ) − H ′ h ( r ) ̺ B (cid:3) u B · n d S x d t ≤ − Z τ Z Ω h ̺ b u · ∂ t U + ̺ b u · ∇ x U · b u + p h ( ̺ )div x U i d x d t + Z τ Z Ω ∂ t p h ( r ) d x d t + Z τ Z Ω h ̺∂ t (cid:18) | U | − H ′ h ( r ) (cid:19) + ̺ b u · ∇ x (cid:18) | U | − H ′ h ( r ) (cid:19) i d x d t + Z τ R Rh [ ̺, v , r, V , u B ]( t )d t + R Eh [ ̺, u ]( τ )where R Rh [ ̺, v , r, V , u B ] = < Π Ch,h , | V | − H ′ h ( r ) > − < Π Mh,h , V > + Z Ω (cid:16) S ( ∇ h u ) : ∇ h (˜ u B − u B ) + p h ( ̺ )div h ( u B − ˜ u B ) + ̺∂ t U · ( b ˜ u B − u B )+ ̺ b u · ∇ h u B · ( b ˜ u B − ˜ u B ) + ̺ b u · ∇ x V · (˜ u B − u B ) + ̺ b u · ∇ h (˜ u B − u B ) · ( V − b v ) (cid:17) d x.
40. Regrouping conveniently several terms in the above expression, we get (cid:20)Z Ω E h (cid:16) ̺, b v (cid:12)(cid:12)(cid:12) r, V (cid:17) d x (cid:21) t = τt =0 + Z τ Z Ω S ( ∇ h u ) : ∇ h ( u − U ) d x d t (10.7)+ Z τ Z Γ out (cid:2) H h ( ̺ ) − H ′ h ( r )( ̺ − r ) − H h ( r ) (cid:3) u B · n d S x d t + Z τ Z Γ in (cid:2) H h ( ̺ B ) − H ′ h ( r )( ̺ B − r ) − H h ( r ) (cid:3) u B · n d S x d t < ∼ − Z τ Z Ω ̺ ( U − b u ) · ∇ x U · ( U − b u ) d x d t − Z τ Z Ω h p h ( ̺ ) − p ′ h ( r )( ̺ − r ) − p h ( r ) i div x U d x d t + Z τ Z Ω ̺r ( U − b u ) · h ∂ t ( r U ) + div x ( r U ⊗ U ) + ∇ x p h ( r ) i d x d t + Z τ Z Ω (cid:16) ̺r ( b u − U ) · U + p ′ h ( r ) (cid:16) − ̺r (cid:17)(cid:17) h ∂ t r + div x ( r U ) i d x d t + Z τ R Rh [ ̺, v , r, V , u B ]( t )d t + R Eh [ ̺, u ]( τ ) ,
4. We shall continue to re-arranging the inequality (10.7) and write (cid:20)Z Ω E h (cid:16) ̺, b v (cid:12)(cid:12)(cid:12) r, V (cid:17) d x (cid:21) t = τt =0 + Z τ Z Ω S ( ∇ h ( u − ˜ U )) : ∇ h ( u − ˜ U ) d x d t (10.8)+ Z τ Z Γ out (cid:2) H h ( ̺ ) − H ′ h ( r )( ̺ − r ) − H h ( r ) (cid:3) u B · n d S x d t + Z τ Z Γ in (cid:2) H h ( ̺ B ) − H ′ h ( r )( ̺ B − r ) − H h ( r ) (cid:3) u B · n d S x d t < ∼ − Z τ Z Ω ̺ ( V − b v ) · ∇ x U · ( V − b v ) d x d t − Z τ Z Ω h p h ( ̺ ) − p ′ h ( r )( ̺ − r ) − p h ( r ) i div x U d x d t + Z τ Z Ω ̺r ( U − b u ) · h ∂ t ( r U ) + div x ( r U ⊗ U ) + ∇ x p ( r ) i d x d t − Z τ Z Ω S ( ∇ x U ) : ∇ h ( u − ˜ U ) d x d t + Z τ Z Ω (cid:16) ̺r ( b u − U ) · U + p ′ h ( r ) (cid:16) − ̺r (cid:17)(cid:17) h ∂ t r + div x ( r U ) i d x d t + Z τ R Rh [ ̺, v , r, V , u B ]( t )d t + R Eh [ ̺, u ]( τ ) , where R Rh [ ̺, v , r, V , u B ] = R Rh [ ̺, v , r, V , u B ] + Z Ω S ( ∇ h u ) : ∇ h ( U − ˜ U ) d x − Z Ω S ( ∇ h ( ˜ U − U )) : ∇ h ( u − ˜ U ) d x + Z Ω ( ̺ b ˜ u B − u B ) · ∇ x U · ( U − b u ) d x + Z Ω ̺ ( b v − V ) · ∇ x V · ( u B − b ˜ u B ) d x + 2 h η Z Ω ̺ ∇ x r · ( U − b u ) d x.
41. The strategy we want to apply recommends to estimate the right hand side of (10.7)–where [ r, U ] is astrong solution of (1.1–1.5)– via the relative energy functional. Seeing the quantities being comparedin the relative energy functional, we shall still rewrite (10.8) as follows: (cid:20)Z Ω E h (cid:16) ̺, b v (cid:12)(cid:12)(cid:12) r, V (cid:17) d x (cid:21) t = τt =0 + Z τ Z Ω S ( ∇ h ( u − ˜ U )) : ∇ h ( u − ˜ U ) d x d t (10.9)+ Z τ Z Γ out (cid:2) H h ( ̺ ) − H ′ h ( r )( ̺ − r ) − H h ( r ) (cid:3) u B · n d S x d t + Z τ Z Γ in (cid:2) H h ( ̺ B ) − H ′ h ( r )( ̺ B − r ) − H h ( r ) (cid:3) u B · n d S x d t < ∼ − Z τ Z Ω ̺ ( V − b v ) · ∇ x U · ( V − b v ) d x d t − Z τ Z Ω h p h ( ̺ ) − p ′ h ( r )( ̺ − r ) − p h ( r ) i div x U d x d t + Z τ Z Ω ̺ ≥ r ( ̺ ) ̺r ( U − b u ) · h ∂ t ( r U ) + div x ( r U ⊗ U ) + ∇ x p ( r ) i d x d t + Z τ Z Ω ̺ ≥ r ( ̺ )(ˆ u − U ) · div x ( S ( ∇ x U )) d x d t + Z τ Z Ω ̺ ≤ r ( ̺ ) ̺r ( ˜ U − u ) · h ∂ t ( r U ) + div x ( r U ⊗ U ) + ∇ x p ( r ) i d x d t + Z τ Z Ω ̺ ≤ r ( ̺ )( u − ˜ U ) · div x ( S ( ∇ x U )) d x d t + Z τ Z Ω (cid:16) ̺r ( b u − U ) · U + p ′ h ( r ) (cid:16) − ̺r (cid:17)(cid:17) h ∂ t r + div x ( r U ) i d x d t + Z τ R Rh [ ̺, v , r, V , u B ]( t )d t + R Eh [ ̺, u ]( τ ) , where R Rh [ ̺, v , r, V , u B ] = R Rh [ ̺, v , r, V , u B ]+ Z Ω ̺ ≥ r ( ̺ )div x ( S ( ∇ x U )) · ( u − ˆ u + U − ˜ U ) d x + Z τ Z Ω ̺ ≤ r ( ̺ ) ̺r ( U − ˜ U + u − b u ) · h ∂ t ( r U ) + div x ( r U ⊗ U ) + ∇ x p ( r ) i d x d t − X K ∈T X σ ∈E ( K ) Z σ n σ,K · S ( ∇ x U ) · ( u − ˜ U )d S x . In the above, we have, among others, integrate by parts in R Ω S ( ∇ x U ) : ∇ h ( u − ˜ U ) d x and then”dispatch” the volume integral conveniently, to the sets { ̺ ≤ r } and { ̺ > r } . We recall that r isdefined in (3.13). 42ince u − ˜ U = v − ˜ V , the latter term in R Rh [ ̺, v , r, V , u B ] is equal to − X K ∈T X σ ∈E ( K ) Z σ n σ,K · S ( ∇ x U ) · ( v − ˜ V )d S x and thus its absolute value is < ∼ √ h . Indeed, for any w ∈ V (Ω; R ) X K ∈T X σ ∈E ( K ) Z σ n σ,K ⊗ ∇ x U ⊗ w d S x = X K ∈T X σ ∈E ( K ) ∩E int Z σ n σ,K ⊗ ( ∇ x U − ( ∇ x U ) σ ) ⊗ ( w − w σ )d S x + X K ∈T X σ ∈E ( K ) ∩E ext Z σ n σ,K ⊗ ∇ x U ⊗ w d S x where we have used the fact that both mean values of w and ∇ x U are continuous over each face σ ∈ E int . Consequently, (cid:12)(cid:12)(cid:12) X K ∈T X σ ∈E ( K ) ∩E int Z σ n σ K ⊗ ∇ x U ⊗ w d S x (cid:12)(cid:12)(cid:12) < ∼ h k∇ x U k L ∞ (Ω) k∇ h w k L (Ω) , (cid:12)(cid:12)(cid:12) X K ∈T X σ ∈E ( K ) ∩E ext Z σ n σ,K ⊗ ∇ x U ⊗ w d S x (cid:12)(cid:12)(cid:12) < ∼ √ h k∇ x U k L ∞ (Ω) k∇ h w k L (Ω) , where we have adapted the reasoning from (8.4) (see also (7.5) and (7.7)).Further, revisiting (10.6), the form of R Rh in (10.8), recalling that R Ω p h ( ̺ )div h ( u B − ˜ u B ) d x = 0, cf.(4.12), using (7.8), (7.22), (7.24) and employing (4.3), (4.6), (4.5) together with estimates proved inLemma 6.2, we deduce | Z τ R Rh [ ̺, v , r, V , u B ]( t )d t | < ∼ Ch α R , R Eh [ ̺, u ]( τ ) < ∼ h, α R = min { α C , α M , η } (10.10)with some number C > r, U ] as indicated in (3.13).We have shown the following result. Lemma 10.1 [ Relative energy inequality for numerical solutions]
Let [ ̺, u ] be a numerical solution of thealgebraic system (3.1–3.4) in the setting (3.6) with the pressure satisfying (1.6–1.8) and the initial andboundary conditions verifying (1.2–1.4). Suppose that the couple [ r, V ] belongs to the class (1.14). Let U = V + u B . Then the relative energy inequality (10.9) holds for all ≤ τ ≤ T . The remainders R Rh , R Eh verify estimates (10.10). Now, we use in the relative energy inequality (10.9) as test functions a strong solution [ r, U ] of problem(1.1–1.5) in the class (1.14).In view of (1.1) we can rewrite (10.9) in the form (cid:20)Z Ω E h (cid:16) ̺, b v (cid:12)(cid:12)(cid:12) r, V (cid:17) d x (cid:21) t = τt =0 + Z τ Z Ω S ( ∇ h ( u − ˜ U )) : ∇ h ( u − ˜ U ) d x d t (10.11)43 Z τ Z Γ out (cid:2) H h ( ̺ ) − H ′ h ( r )( ̺ − r ) − H h ( r ) (cid:3) u B · n d S x d t < ∼ − Z τ Z Ω ̺ ( V − b v ) · ∇ x U · ( V − b v ) d x d t − Z τ Z Ω h p h ( ̺ ) − p ′ h ( r )( ̺ − r ) − p h ( r ) i div x U d x d t + Z τ Z Ω ̺ ≤ r ( ̺ ) ̺ − rr ( ˜ U − u ) · div x S ( ∇ x U ) d x d t + Z τ Z Ω ̺ ≥ r ( ̺ ) ̺ − rr ( U − b u ) · div x S ( ∇ x U ) d x d t + h α R , where we have used the inequality (cid:12)(cid:12)(cid:12) Z τ Z Γ in (cid:2) H h ( ̺ B ) − H ′ h ( r )( ̺ B − r ) − H h ( r ) (cid:3) u B · n d S x d t (cid:12)(cid:12)(cid:12) < ∼ h (recall that ̺ B = b r B ).At this stage, it is convenient to recall a simple but in our context important algebraic lemma: Lemma 10.2
Let < a < b < ∞ and let p satisfiy (1.6–1.8). Then there exists a number c = c ( a, b ) > such that for all ̺ ∈ [0 , ∞ ) and r ∈ [ a, b ] , E ( ̺ | r ) ≥ c ( a, b ) (cid:16) O res ( ̺ ) + ̺ O res ( ̺ ) + 1 O res ( ̺ ) p ( ̺ ) + ( ̺ − r ) O ess ( ̺ ) (cid:17) , (10.12) where E is defined in (1.9) and O ess = [ a/ , b ] , O res = [0 , ∞ ) \ O ess . With Lemma 10.2 it is rudimentary to see that (cid:12)(cid:12)(cid:12) Z τ Z Ω ̺ ( V − b v ) · ∇ x U · ( V − b v ) d x d t (cid:12)(cid:12)(cid:12) < ∼ C Z τ E ( ̺, b v | r, V )d t and, under assumptions (1.6–1.8) on the pressure, also (cid:12)(cid:12)(cid:12) Z τ Z Ω h p ( ̺ ) − p ′ ( r )( ̺ − r ) − p ( r ) i div x U d x d t (cid:12)(cid:12)(cid:12) < ∼ C Z τ E ( ̺, b v | r, V )d t, We observe that, under (1.6–1.8), also (cid:12)(cid:12)(cid:12) Z τ Z ̺ ≥ r ̺ − rr ( U − ˆ u ) · div x S ( ∇ x U )d x d t (cid:12)(cid:12)(cid:12) < ∼ C (cid:16) Z τ E ( ̺, b v | r, V )d t + h (cid:17) while (cid:12)(cid:12)(cid:12) Z τ Z ̺ ≤ r ̺ − rr ( ˜ U − u ) · div x S ( ∇ x U )d x d t (cid:12)(cid:12)(cid:12) < ∼ δ Z τ k ˜ U − u k L (Ω) d t + c ( δ ) C Z τ E ( ̺, b v | r, V )d t In the last four estimates δ > c = c ( δ ) > C > r is defined in (3.13).Finally recalling the inequality, cf. (4.11), k ˜ U − u k L (Ω) < ∼ k∇ h ( ˜ U − u ) k L (Ω) , (cid:20)Z Ω E h (cid:16) ̺, b v (cid:12)(cid:12)(cid:12) r, V (cid:17) d x (cid:21) t = τt =0 + Z τ (cid:16) k ˜ U − u k L (Ω) + k∇ h ( ˜ U − u ) k L (Ω) (cid:17) d t < ∼ C (cid:16) h α R + Z τ E h ( ̺, b v | r, V )d t (cid:17) . Applying to (10.11) the Gronwall argument, we get the following lemma
Lemma 10.3 [Error estimates for the numerical solution]
Let assumptions of Lemma 10.1 be satisfied. Let [ ̺ h , u h ] be a numerical solution of the algebraic system (3.1–3.4) in the setting (3.6) emanating from theinitial and boundary value data (1.2–1.4). Suppose that the couple [ r, U ] , U = V + u B belongs to the class(1.14) and represents a strong solution of the problem (1.1–1.5). Then (cid:20)Z Ω E h (cid:16) ̺ h , b v h (cid:12)(cid:12)(cid:12) r, V (cid:17) d x (cid:21) t = τt =0 + Z τ (cid:16) k ˜ U − u h k L (Ω) + k∇ h ( ˜ U − u h ) k L (Ω) (cid:17) d t < ∼ Ch α R , where α R = min { α C , α M , η } > with C > dependent on the strong solution in the way indicated in (10.7). This finishes the proof of Theorem 3.2.
Remark 10.1
In view of estimates evoked in Remark 6.2, and due to Remarks 7.1, 7.2 and (7.24), Lemma10.3 remains true with α R = min { α C , α M } > also for the numerical solutions of the scheme (3.1–3.4)with ˜ κ = 0 and p in the class (1.6–1.8) in the situations described by conditions (3.15),(3.16) resp. (3.15),(3.17), (3.18), resp. (3.15), (3.17), (3.19). We let the details for the interested reader.
11 Concluding remarks
In this section, we wish to mention some open problems related to the numerical approximations of theNavier-Stokes equations with non-homogenous boundary data.1. In this paper, we consider only the case when the computational domain Ω h = ∪ K ∈T h K h coincideswith the physical domain Ω. If this is not the case, one must work on so called unfitted meshes (seee.g. Babuˇska [2]), i.e. Ω h = Ω. For the no-slip boundary conditions, a solution has been proposed in[19], [15] but for the non-homogenous boundary data, this problem is clearly more involved. Namely,the quality of approximation of Ω by Ω h and of ∂ Ω by ∂ Ω h will play a preponderant role, and remainsto be determined.2. Karper’s scheme [50] is a collocalized scheme (the density and velocity are discretized on the same”primal” mesh). Quite often, the staggered schemes (the density resp. the pressure are discretized ona primal mesh while the velocity components on a ”bi-dual” meshes) –as e.g. the Marker and Cell(MAC) finite difference schemes, [41], [37], [39], [31], [35] or staggered discretizations combined withRannacher-Turek finite elements or Crouzeix-Raviart finite elements [43]– are computationally moreefficient. It would be certainly of interest to extend the convergence results to dissipative solutions tothese types of schemes. 45. Since the proof of convergence to dissipative solutions is relatively weakly scheme-structure dependentit would be of great interest to determine as weak as possible universal criteria imposed on the numericalscheme in order to obtain the convergence to the dissipative solutions. The universal description ofseveral staggered schemes in one unique formalism provided in [43] may be a starting point to approachthis task.4. Gallouet et al. established in [38] the error estimates for the Karper’s scheme and in [39] for the MACscheme under assumption γ > / t and the sizeof the space discretization h by using a different approach, which consists in mimicking the proof ofthe weak-strong uniqueness known from the continuous case (see [18]). For ”large” values of γ , theseestimates seem to be optimal judging from what can be obtained for the sole continuity equation withthe same regularity of the transporting velocity field, see [54] for the discussion about the optimality. Itis certainly of interest to try to prove similar results with the general non homogenous boundary data.The weak strong uniqueness principle established in this boundary value setting in the continuous casein [52] and later in [1] encourages the attempts in this direction.5. Weak solutions to the problem (1.1–1.5) are known to exist for the adiabatic coefficients γ > /
2, see[25] for the no-slip case and [6], [7], [52] for the general inflow-outflow boundary conditions. The proofsuse, among others, very specific mathematical tools closely related to the structure of the system –as. e.g. compensated compactness and various commutator lemmas– the tools whose numericalcounterpart is usually not available. The consistency formulation of the balance laws allows to bring thenumerical solution close to a weak solution modulo a remainder, which in the case of the convergenceto weak solutions must have a ”better quality” than in the proofs of the convergence to the dissipativesolutions. T. Karper provided in 2013 in [50] a convergence proof of the Karper’s scheme to a weaksolution under assumption γ >
References [1] A. Abbatiello, E. Feireisl, and A. Novotn´y. Generalized solutions to models of compressible viscousfluids. Arxiv preprint No. 1912.12896, 2019.[2] I. Babuˇska. The finite element method with penalty
Math. Comput.
27: 221-228, 1973.[3] J. W. Barrett, E. S¨uli. Existence of global weak solutions to compressible isentropic finitely extensiblenonlinear bead-spring chainmodels for dilute polymers.
Math. Models Methods Appl. Sci.
16 (2016)DOI: 10.1142/S0218202516500093[4] M. Bessemoulin-Chatard, C. Chainais-Hillairet, F. Filbet On discrete functional inequalities for somefinite volume schemes.
IMA Journal of Numerical Analysis , 2014465] F. Brezzi, M. Fortin.
Mixed and hybrid finite element methods . Springer Series in ComputationalMathematics. Vol. 15, Springer-Verlag, New York, 1991.[6] T. Chang, B. J. Jin, and A. Novotn´y. Compressible Navier-Stokes system with general inflow-outflowboundary data.
SIAM J. Math. Anal. , (2):1238–1278, 2019.[7] H.J. Choe, A. Novotny, M. Yang Compressible Navier-Stokes system with general inflow-outflow bound-ary data on piecewise regular domains ZAMM Z. Angew. Math. Mech.
Rev. Francaise Automat. In-format. Recherche Op´erationnelle S´er. Rouge
7: 3375, 1973[9] V. Dolejˇsi, E. Feistauer, J. Felcman, A. Klikov´a. Error estimates for barycentric finite volumes combinedwith nonconforming finite elements applied to nolinear convection diffusion problems.
Appl. Maths.
Partial Differential Equations.
Graduate Studies in Mathematics, vol. 19, AMS Providence,1998.[11] R. Eymard, T. Gallou¨et, and R. Herbin. Finite volume methods.
Handbook of numerical analysis ,7:713–1018. North-Holland, Amsterdam, 2000.[12] R. Eymard, T. Gallou¨et, R. Herbin J-C. Latch´e. A convergent finite element-finite volume scheme forthe compressible Stokes problem. II. The isentropic case.
Math. Comp. , 270(79):649–675, 2010.[13] E. Feireisl.
Dynamics of viscous compressible fluids , volume 26 of
Oxford Lecture Series in Mathematicsand its Applications . Oxford University Press, Oxford, 2004.[14] E. Feireisl, P. Gwiazda, A. Swierczewska-Gwiazda, E Wiedemann Dissipative measurevalued solutions to the compressible Navier-Stokes system.
Calc. Var.
55, 141 (2016). https://doi.org/10.1007/s00526-016-1089-1 [15] E. Feireisl, T. Karper, M. Mich´alek. Convergence of a numerical method for the com-pressible Navier-stokes system on general domains.
Numer. Math. https://doi.org/10.1007/s00211-015-0786-6 [16] E. Feireisl, T. Karper, M. Pokorn´y.
Mathematical theory of compressible viscous fluids: Analysis andNumerics.
Springer. Heidelberg. 2016[17] E. Feireisl, T. Karper, A. Novotn´y. A convergent numerical method for the Navier-Stokes-Fouriersystem
IMA J.Numer. Anal.
36: 1477-1535, 2016[18] E. Feireisl, B J. Jin, and A. Novotn´y. Relative entropies, suitable weak solutions, and weak-stronguniqueness for the compressible Navier–Stokes system.
J. Math. Fluid. Mech. , 14(4):717–730, 2012.[19] E. Feireisl, R. Hoˇsek, D. Maltese, and A. Novotn´y. Error estimates for a numerical method for thecompressible Navier-Stokes system on sufficiently smooth domains.
ESAIM : Mathematical Modelingand Numerical Analysis , 51(1):279–319, 2017. 4720] E. Feireisl, Y. Lu, E. S¨uli Dissipative weak solutions to compressible NavierStokesFokkerPlanck systemswith variable viscosity coefficients.
Journal of Mathematical Analysis and Applications , 443 (1):322-351,2016.[21] E. Feireisl, M. Luk´aˇcov´a-Medvidov´a. Convergence of a mixed finite element finite volume scheme forthe isentropic Navier-Stokes system via dissipative measure-valued solutions.
Found Comput Math https://doi.org/10.1007/s10208-017-9351-2 [22] E. Feireisl, M. Luk´aˇcov´a-Medvidov´a, H. Mizerov´a A finite volume scheme for the Euler system inspiredby he two velocities approach. Arxiv 1805.05072, 2019[23] E. Feireisl, M. Luk´aˇcov´a-Medvidov´a, H. Mizerov´a, B. She. Convergence of a finite volume scheme forthe compressible Navier-Stokes syastem Arxiv 1811.02866, 2019[24] E. Feireisl, M. Luk´aˇcov´a-Medvid´ov´a, and S. Neˇcasov´a, A. Novotn´y, and B. She. Asymptotic preservingerror estimates for numerical solutions of compressible Navier-Stokes equations in the low Mach numberregime.
Multiscale Model. Simul. , 16(1), 150–183, 2017[25] E. Feireisl, A. Novotn´y and H. Petzeltov´a,
On the existence of globally defined weak solutions to theNavier-Stokes equations of compressible isentropic fluids , J. Math. Fluid Mech., 3(2001), 358–392.[26] M. Feistauer.
Mathematical methods in fluid dynamics.
Pitman Monographs and Surveys in Pure andApplied Mathematics Series 67, Longman, Harlow, 1993[27] M. Feistauer, J. Felcman, M. Luk´aˇcov´a-Medvidov´a. Combined finite element-finite volume solution ofcompressible flow.
J. Comput. Appl. Math.
63 : 179–199, 1995[28] M. Feistauer, J. Felcman, M. Luk´aˇcov´a-Medvidov´a, G. Warnecke. Error estimates of a combined finiteelement-finite volume solution for non linear convection-diffusion problems.
SIAM J. Numer. Anal.
Mathematical and computational methods for compressibleflow.
Clarendon Press. Oxford, 2003.[30] M. Feistauer, J. Cesenek, V. Kuˇcera. Discontinuous Galerkin method– a robust solver for compressibleflow.
Recent developments in the numerics of conservation laws. Notes Numer. Fluid. Multidiscip. Des.
IMAJ. Numer. Anal. , 33(3):922–951, 2013.[32] U. S. Fjordholm, S. Mishra and E. Tadmor. Arbitrarily high order accurate entropy stable essentiallynon-oscillatory schemes for systems of conservation laws.
SIAM J. Numer. Anal.
Found. Comput. Math.
ActaNumer.
M2AN Math. Model. Numer. Anal.
Math. Comp. , 267(78):1333–1352, 2009.[37] T. Gallou¨et, R. Herbin, J-C Latch´e, and K Mallem. Convergence of the MAC scheme for the in-compressible Navier-Stokes equations.
Foundations of Computational Mathematics , 18(1): 249–289,2018.[38] T. Gallou¨et, R. Herbin, D. Maltese, and A. Novotny. Error estimates for a numerical approximation tothe compressible barotropic Navier–Stokes equations.
IMA J. Numer. Anal. , 36(2):543-592, 2016.[39] T. Gallou¨et, D. Maltese, and A. Novotn´y. Error estimates for the implicit MAC scheme for the com-pressible Navier-Stokes equations.
Numer. Math. , 141(2):495–567, 2019.[40] V. Girinon. Navier-Stokes equations with nonhomogeneous boundary conditions in a bounded three-dimensional domain.
J. Math. Fluid Mech. , 13: 309–339, 2011[41] F.H. Harlow and A.A. Amsden. A numerical fluid dynamics calculation method for all flow speeds.
J.Comp. Phys. , 8:197–213, 1971.[42] R. Herbin, W. Kheriji, J-C. Latch´e On some implicit and semi-implicit staggered schemes for theshallow water and Euler equations. ESAIM Math. Model. Numer. Anal. 48(6): 1807–1857, 2014.[43] R. Herbin, J-C. Latch´e, S. Minjeaud and N. Therme. Conservativity and weak consistency of a class ofstaggered finite volume methods for the Euler equations. https://arxiv.org/pdf/1910.03998.pdf ,2019[44] R. Hosek, B. She Stability and consistency of a finite difference scheme for compressible viscousisentropic flow in multi-dimension.
J. Numer. Math.
Kragujevac J. Math , 30(1):263–275, 2007.[46] V. Jovanovi´c and C. Rohde. Finite-volume schemes for friedrichs systems in multiple space dimen-sions: A priori and a posteriori error estimates.
Numerical Methods for Partial Differential Equations ,21(1):104–131, 2005.[47] K. H. Karlsen, T. Karper. A convergent nonconforming finite element method for compressible Stokesflow.
SIAM J. Numer. Anal.
Math. Comp.
IMA J. Numer. Anal.
Numer.Math. , 125(3) : 441–510, 2013. 4951] D. Kr¨oner
Numerical schemes for conservation laws.
Willey-Teubner, 1997[52] Y. Kwon, A. Novotny. Dissipative solutions to compressible Navier-Stokes equations with generalinflow-outflow data: existence, stability and weak-strong uniqueness. arXiv:1905.02667, 2019[53] P-L. Lions.
Mathematical topics in fluid mechanics. Vol. 2 , volume 10 of
Oxford Lecture Series inMathematics and its Applications . Oxford University Press, New York, 1998.[54] D. Maltese, A. Novotny. Implicit MAC scheme for compressible Navier-Stokes equations: Low Machnumber asymptotics. IMA J. Numer. Anal. 2020, accepted[55] A. Novotny. Weak solutions for a bi-fluid model for a mixture of two compressible non interactingfluids.
Sci. China Math. , 2019, DOI10.1007/s11425-019-9552-1.[56] A. Novotn´y and I. Straskraba. Introduction to the mathematical theory of compressible flow. volume 27of
Oxford Lecture Series in Mathematics and its Applications . Oxford University Press, New York, 2004.[57] R. Rannacher and S. Turek. Simple nonconforming quadrilateral Stokes element.
Numerical Methodsfor Partial Differential Equations , 8(2):97–111, 1992.[58] E. Tadmor Perfect derivatives, conservative differences and entropy stable computation of hyperbolicconservation laws.
Dicr. Cont. Systems
Acta Numer.
12: 451–512, 2003[60] E. Tadmor, W. Zhong Entropy stable approximations of Navier-Stokes equations with no artificialnumerical viscosity.
J. Hyperbolic Diff. Eq.
3: 529–559, 2006[61] A. Valli and W. M. Zajaczkowski. Navier-Stokes equations for compressible fluids: global existence andqualitative properties of the solutions in the general case.