On the influence of the nonlinear term in the numerical approximation of Incompressible Flows by means of proper orthogonal decomposition methods

Available online at www.sciencedirect.com ScienceDirect Comput. Methods Appl. Mech. Engrg. 405 (2023) 115866 www.elsevier.com/locate/cma On the influence of the nonlinear term in the numerical approximation of Incompressible Flows by means of proper orthogonal decomposition methods Bosco García-Archillaa,1, Julia Novob,∗,2, Samuele Rubinoc,3 aDepartamento de Matemática Aplicada II, Universidad de Sevilla, Sevilla, Spain bDepartamento de Matemáticas, Universidad Autónoma de Madrid, Spain cDepartment EDAN & IMUS, Universidad de Sevilla, Spain Received 29 March 2022; received in revised form 19 December 2022; accepted 20 December 2022 Available online 30 December 2022 Abstract We consider proper orthogonal decomposition (POD) methods to approximate the incompressible Navier–Stokes equations. We study the case in which one discretization for the nonlinear term is used in the snapshots (that are computed with a full order method (FOM)) and a different discretization of the nonlinear term is applied in the POD method. We prove that an additional error term appears in this case, compared with the case in which the same discretization of the nonlinear term is applied for both the FOM and the POD methods. However, the added term has the same size as the error coming from the FOM so that the rate of convergence of the POD method is barely affected. We analyze the case in which we add grad–div stabilization to both the FOM and the POD methods because it allows to get error bounds with constants independent of inverse powers of the viscosity. We also study the case in which no stabilization is added. Some numerical experiments support the theoretical analysis. © 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). MSC: 35Q30; 65M12; 65M15; 65M20; 65M60; 65M70; 76B75 Keywords: Navier–Stokes equations; Proper orthogonal decomposition; Nonlinear term discretization; Grad–div stabilization 1. Introduction The computational cost of direct numerical simulations can be reduced by using reduced order models. The proper orthogonal decomposition (POD) method is based on a reduced basis that is computed using snapshots coming from a full order method (FOM). In this paper, we study the numerical approximation of incompressible flows with POD methods. ∗Corresponding author. E-mail addresses: [email protected] (B. Garc´ ıa-Archilla), julia.nov[email protected] (J. Novo), [email protected] (S. Rubino). 1Research is supported by Spanish MCINYU under grants PGC2018-096265-B-I00 and PID2019-104141GB-I00. 2Research is supported by Spanish MINECO under grants PID2019-104141GB-I00 and VA169P20. 3Research is supported by Spanish MCINYU under grant RTI2018-093521-B-C31 and by the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie Actions, grant agreement 872442-ARIA. https://doi.org/10.1016/j.cma.2022.115866 0045-7825/© 2022 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http: //creativecommons.org/licenses/by-nc-nd/4.0/). B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 We consider the Navier–Stokes equations ∂tu−ν∆u+(u·∇)u+∇p=fin (0,T]×Ω, ∇ ·u=0 in (0,T]×Ω,(1) in a bounded domain Ω⊂Rd,d∈ {2,3}with initial condition u(0) =u0. In (1),uis the velocity field, pthe kinematic pressure, ν > 0 the kinematic viscosity coefficient, and frepresents the accelerations due to external body forces acting on the fluid. The Navier–Stokes equations (1) must be complemented with boundary conditions. For simplicity, we only consider homogeneous Dirichlet boundary conditions u=0on ∂Ω. In [1] we can find a definition of FOM–ROM consistency. A ROM is said FOM consistent if it uses the same computational model and numerical discretization as the FOM. As stated in [1]: the theoretical investigation of the FOM–ROM consistency is scarce and investigating the FOM–ROM consistency both computationally and theoretically is an important research direction in reduced order modeling. As in [1], in the present paper we study the consistency related to the nonlinear term. In practical simulations one can apply some given software to compute the snapshots. It could then be the case that a different discretization is used for the discretization of the nonlinear term in the FOM method and the POD method (the last one being typically implemented by means of a hand-made code instead of an existing one). Our aim in this paper is to study the influence of the use of different discretizations for the nonlinear terms on the final error bounds of the POD method. In [2] POD stabilized methods for the Navier–Stokes equations were considered and analyzed for a case in which the snapshots are based on a non inf–sup stable method and a case in which the snapshots are based on an inf–sup stable method. In the present paper, we consider the second case. Our snapshots are based on an inf–sup stable method and, as in [2], we add grad–div stabilization to both the FOM and the POD methods. We analyze the case in which different discretizations are applied for the nonlinear term. Adding grad–div stabilization we are able to prove error bounds for the method with constants that do not depend explicitly on inverse powers of the viscosity, although, as usual, may depend on it through the norms of the theoretical solution. In the last section, devoted to numerical experiments, we also study the use of different values for the grad–div parameter in the FOM and ROM models. We conclude that it is advisable to use the same values for the grad–div parameter in both models. In the recent paper [1] we have found an study of the consistency of the nonlinear discretization in FOM and POD methods. In this reference no stabilization is included neither in the FOM nor in the POD method. The authors of [1] conclude that the use of different discretizations in the nonlinear term yields additional terms that prevent the POD method from recovering the FOM accuracy. In the present paper, we prove, in agreement with the results in [1], that there are some additional terms in the error bound of the final POD method coming from the use of different discretizations in the nonlinear terms. We prove that the additional terms have the size of the L2error in the velocity and the L2error of the divergence of the velocity of the FOM method. This holds both for the methods with grad–div stabilization and also for the plain methods considered in [1]. We then conclude that the extra error term coming from using different discretization for the nonlinear terms in the FOM and POD methods leads in general to a less accurate method although in some examples the extra error is still small and the non consistent FOM–ROM method could be used in practice. Following [1], we analyzed a concrete case in which the so called divergence or skew-symmetric form of the nonlinear term (commonly used in practice) is used for the FOM method and the EMAC form, see [3], is applied for the POD method. The error analysis for any other combination of discretizations of the nonlinear terms with analogous properties could be carried out in a similar way. The EMAC formulation was designed to conserve energy, momentum and angular momentum. Considering a continuous in time method it is easy to see, as we prove in Section 4, that the effect of grad–div stabilization in a method with EMAC form for the nonlinear term is a lost of kinetic energy. The method maintains however the conservation of the momentum and angular momentum. In practice, as stated in [4], a fully version of a method with EMAC form for the nonlinear term is applied. In [4] it is proved that using Newton method for the nonlinear term is the option from which more quantities are still conserved, although the conservation of the kinetic energy is lost. In this sense, adding grad–div stabilization incurs in the same phenomena so that it seems not to be a hard problem since either kinetic energy or any other properties will be lost with the EMAC formulation in the fully discrete case. The outline of the paper is as follows. In Section 2we introduce some notation. In Section 3we state some preliminaries concerning the POD method. Section 4is devoted to the error analysis in which one considers different discretizations for the nonlinear term in the FOM and POD methods. In Section 5we present some numerical experiments. We end the paper with some conclusions. 2 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 2. Preliminaries and notation The following Sobolev embeddings [5] will be used in the analysis: For q∈[1,∞), there exists a constant C=C(Ω,q) such that ∥v∥Lq′≤C∥v∥Ws,q,1 q′≥1 q−s d>0,q<∞, v ∈Ws,q(Ω)d.(2) The following inequality can be found in [6, Remark 3.35] ∥∇ ·v∥0≤ ∥∇v∥0,v∈H1 0(Ω)d,(3) where, here and in the sequel, we use the notation ∥·∥jfor ∥·∥Wj,2= ∥· ∥Hj. Let us denote by Q=L2 0(Ω)= {q∈L2(Ω)|(q,1) =0}. Let Th=(τh j, φh j)j∈Jh,h>0 be a family of partitions of Ω, where hdenotes the maximum diameter of the elements τh j∈Th, and φh jare the mappings from the reference simplex τ0onto τh j. We shall assume that the partitions are shape-regular and quasi-uniform. We define the following finite element spaces Yl h={vh∈C0(Ω)|vh|K∈Pl(K),∀K∈Th},l≥1, Yl h=(Yl h)d,Xl h=Yl h∩H1 0(Ω)d, Ql h=Yl h∩L2 0(Ω). Vh,l=Xl h∩{χh∈H1 0(Ω)d|(qh,∇ ·χh)=0∀qh∈Ql−1 h},l≥2.(4) If the family of meshes is quasi-uniform then the following inverse inequality holds for each vh∈Yl h, see e.g., [7, Theorem 3.2.6], ∥vh∥Wm,p(K)≤cinvhn−m−d(1 q−1 p) K∥vh∥Wn,q(K),(5) where 0 ≤n≤m≤1, 1 ≤q≤p≤ ∞, and hKis the diameter of K∈Th. Let l≥2, we consider the MFE pair known as Hood–Taylor elements [8,9] (Xl h,Ql−1 h). For these elements a uniform inf–sup condition is satisfied (see [8]), that is, there exists a constant βis >0 independent of the mesh size hsuch that inf qh∈Ql−1 h sup vh∈Xl h (qh,∇ ·vh) ∥vh∥1∥qh∥L2/R≥βis.(6) We will denote by PQhthe L2orthogonal projection onto Ql−1 h. The following bound holds ∥q−PQhq∥0≤Chl∥q∥l,∀q∈Hl(Ω).(7) As a direct method, as in [2, Section 5], we consider a Galerkin method with grad–div stabilization and for simplicity in the analysis we consider the implicit Euler method in time. Let us fix T>0 and M>0 and take ∆t=T/M. The method reads as follows: given u0 h≈u0find (un h,pn h)∈Xl h×Ql−1 hfor n≥1 such that (un h−un−1 h ∆t,vh)+ν(∇un h,∇vh)+bh(un h,un h,vh)−(pn h,∇ ·vh) +µ(∇ ·un h,∇ ·vh)=(fn,vh)∀vh∈Xl h,(8) (∇ ·un h,qh)=0∀qh∈Ql−1 h, where µis the positive grad–div stabilization parameter. For the discretization of the nonlinear term we consider the following form, the so called divergence or skewsymmetric form, bh(u,v,w)=((u·∇v),w)+1 2((∇ ·u)v,w).(9) 3 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 It is well-known that considering the discrete divergence-free space Vh,lwe can remove the pressure from (8) since un h∈Vh,lsatisfies for n≥1 (un h−un−1 h ∆t,vh)+ν(∇un h,∇vh)+bh(un h,un h,vh) +µ(∇ ·un h,∇ ·vh)=(fn,vh),∀vh∈Vh,l.(10) For this method the following bound holds, see [10] ∥un−un h∥0+⎛ ⎝µ M ∑ j=1 ∆t∥∇ ·(un−un h)∥2 0⎞ ⎠ 1/2 ≤C(u,p,l+1) (hl+∆t),(11) for 1 ≤n≤M, where the constant C(u,p,l+1) depends on ∥u∥L∞(Hl+1),(∫T 0∥ut∥2 l)1/2 ,(∫T 0∥utt ∥2 0)1/2 and ∥p∥L∞(Hl)but does not depend explicitly on inverse powers of ν. For the plain Galerkin method the following bound holds where the error constants depend explicitly on inverse powers of ν ∥un−un h∥0+h∥∇(un−un h)∥0≤C(u,p, ν−1,l+1) (hl+1+∆t).(12) 3. Proper orthogonal decomposition We will consider a proper orthogonal decomposition (POD) method. As for the FOM we fix T>0 and M>0 and take ∆t=T/M. We consider the following space U= ⟨u1 h,...,uM h⟩, where uj h=uh(·,tj). Let dvbe the dimension of the space U. Let Kvbe the correlation matrix corresponding to the snapshots Kv=((kv i,j)) ∈RM×M, where kv i,j=1 M(ui h,uj h), and (·,·) is the inner product in L2(Ω)d. Following [11] we denote by λ1≥λ2, . . . ≥λdv>0 the positive eigenvalues of Kvand by v1,...,vdv∈RMthe associated eigenvectors. Then, the (orthonormal) POD bases are given by ϕk=1 √M 1 √λk M ∑ j=1 vj kuh(·,tj),(13) where vj kis the jth component of the eigenvector vkand the following error formulas hold, see [11, Proposition 1] 1 M M ∑ j=0∥uj h− r ∑ k=1 (uj h,ϕk)ϕk∥2 0≤ dv ∑ k=r+1 λk.(14) We will denote by Svthe stiffness matrix for the POD basis: Sv=((sv i,j)) ∈Rdv×dv,sv i,j=(∇ϕi,∇ϕj). In that case, for any v∈U, the following inverse inequality holds, see [11, Lemma 2] ∥∇v∥0≤√∥Sv∥2∥v∥0.(15) In the sequel we will denote by Ur= ⟨ϕ1,ϕ2,...,ϕr⟩, and by Pv r, the L2-orthogonal projection onto Ur. 4 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 3.1. A priori bounds for the orthogonal projection onto Ur In [2, Section 3.1] some a priori bounds for the FOM and for the orthogonal projection were obtained assuming for the solution of (1) the same regularity needed to prove (11). For 0 ≤j≤Mthe following bound holds (see [2, (31), (35)]) ∥uj h∥∞≤Cu,inf (16) ∥∇uj h∥L2d/(d−1) ≤Cu,ld.(17) For 1 ≤j≤Mthe following bounds hold (see [2, (39), (40), (41)]) ∥Pv ruj h∥∞≤Cinf,(18) ∥∇Pv ruj h∥∞≤C1,inf,(19) ∥∇Pv ruj∥L2d/(d−1) ≤Cld.(20) 4. The POD method We now consider the grad–div POD model. As for the full order model, for simplicity in the error analysis, we also use the implicit Euler method as time integrator. Taking an initial approximation u0 r∈Ur, for n≥1, find un r∈Ursuch that (un r−un−1 r ∆t,ϕ)+ν(∇un r,∇ϕ)+bpod (un r,un r,ϕ)+µ(∇ ·un r,∇ ·ϕ) =(fn,ϕ),∀ϕ∈Ur,(21) where we observe that Ur⊂Vh,lso that there is no pressure approximation in (21). For the discretization of the nonlinear term we consider the EMAC form bpod (u,v,w)=(2D(u)v,w)+((∇ ·u)v,w),(22) where the deformation tensor is D(u)=1 2(∇u+(∇u)T), and the second term in (22) is added to ensure the property bpod (v,v,v)=0.(23) The EMAC form is based on the identity (u·∇)u=2D(u)u−1 2∇|u|2.(24) From (24) it can be observed that the EMAC form of the nonlinear term implies a modification of the pressure, i.e., in a velocity–pressure formulation of a method with the EMAC form for the nonlinear term the pressure approximation converges to p−1 2|u|2instead of to the original pressure pin (1). Considering for simplicity the continuous in time case, i.e., the method: (ur,t,ϕ)+ν(∇ur,∇ϕ)+bpod (ur,ur,ϕ)+µ(∇ ·ur,∇ ·ϕ)=(f,ϕ), taking ϕ=urand assuming as in [3,4], f=0 and ν=0 one gets applying (23) d dt ∥ur∥2 0= −µ∥∇ ·ur∥2 0.(25) This means that the grad–div stabilization produces a lost in the kinetic energy. As stated in the introduction, considering a fully discrete method with Newton discretization for the nonlinear term one has also the effect of a lost of kinetic energy, see [4]. On the other hand, it is easy to prove that the grad–div term has no negative effect in the conservation of linear and angular momentum so that both are still conserved. To prove this, one can argue as in [12, Theorem 3.3] and observe that the test functions used to achieve the conservation of linear and angular momentum have divergence zero so that the grad–div term does not affect the proof. 5 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 4.1. Error analysis of the method Denoting by ηn h=Pv run h−un h, it is easy to get (Pv run h−Pv run−1 h ∆t,ϕ)+ν(∇Pv run h,∇ϕ)+bh(Pv run h,Pv run h,ϕ) +µ(∇ · Pv run h,∇ ·ϕ) =(fn,ϕ)+ν(∇ηn h,∇ϕ)+µ(∇ ·ηn h,∇ ·ϕ) (26) +bh(Pv run h,Pv run h,ϕ)−bh(un h,un h,ϕ),∀ϕ∈Ur. Subtracting (26) from (21) and denoting by en r=un r−Pv run h we get ∀ϕ∈Ur (en r−en−1 r ∆t,ϕ)+ν(∇en r,∇ϕ)+µ(∇ ·en r,∇ ·ϕ) (27) =(bh(un h,un h,ϕ)−bpod (un r,un r,ϕ))−ν(∇ηn h,∇ϕ)−µ(∇ ·ηn h,∇ ·ϕ). Taking ϕ=en rwe obtain 1 2∆t(∥en r∥2 0−∥en−1 r∥2 0)+ν∥∇en r∥2 0+µ∥∇ ·en r∥2 0(28) ≤(bh(un h,un h,en r)−bpod (un r,un r,en r))−ν(∇ηn h,∇en r)−µ(∇ ·ηn h,∇ ·en r) =I+I I +III. We will bound the first term on the right-hand side above. Arguing as in [1] we get |I| ≤ ⏐⏐bpod (un h,un h,en r)−bpod (un r,un r,en r)⏐⏐ +⏐⏐bh(un h,un h,en r)−bpod (un h,un h,en r)⏐⏐.(29) To bound the first term on the right-hand side of (29) we argue as in [13] (see also [14]). We first observe that bpod (un h,un h,en r)−bpod (un r,un r,en r)=bpod (−ηn h,un h,en r) +bpod (Pv run h,−ηn h,en r)−bpod (Pv run h,en r,en r) −bpod (en r,Pv run h,en r)−bpod (en r,en r,en r).(30) The last term on the right-hand side of (30) vanishes due to the property (23). The expression of the fourth one is bpod (en r,Pv run h,en r)=((Pv run h·∇)en r,en r)+((en r·∇)en r,Pv run h)+((∇ ·en r)en r,Pv run h). Integration by parts reveals that the first term on the right-hand side above equals −1 2((∇ · Pv run h)en r,en r), whilst the other two equal −((en r·∇)Pv run h,en r), so that bpod (en r,Pv run h,en r)= −1 2((∇ · Pv run h)en r,en r)−((en r·∇)Pv run h,en r). Since the third term on the right-hand side of (30) is bpod (Pv run h,en r,en r)=2((en r·∇)Pv run h,en r)+((∇ · Pv run h)en r,en r), 6 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 we see that the last three terms on the right-hand side of (30) add to −((en r·∇)Pv run h,en r)−1 2((∇ · Pv run h)en r,en r), so that (30) can be written as bpod (un h,un h,en r)−bpod (un r,un r,en r)=bpod (−ηn h,un h,en r) (31) +bpod (Pv run h,−ηn h,en r)−1 2((∇ · Pv run h)en r,en r)−(D(Pv run h)en r,en r). For the first term on the right-hand side of (31), applying (16), we obtain bpod (−ηn h,un h,en r)≤C∥ηn h∥1∥un h∥∞∥en r∥0≤CCu,inf∥ηn h∥1∥en r∥0.(32) Applying H¨ older’s inequality, (19) and Sobolev embedding (2) we get for the second term on the right-hand side of (31) bpod (Pv run h,−ηn h,en r)≤2∥D(Pv run h)∥L2d/(d−1) ∥ηn h∥L2d∥en r∥0 +∥∇ · Pv run h∥L2d/(d−1) ∥ηn h∥L2d∥en r∥0 ≤C∥∇Pv run h∥L2d/(d−1) ∥ηn h∥L2d∥en r∥0 ≤CCld∥ηn h∥1∥en r∥0.(33) Finally, for the last two terms on the right-hand side of (31), applying (19), we get −1 2((∇ · Pv run h)en r,en r)−(D(Pv run h)en r,en r)≤CC1,inf∥en r∥2 0.(34) From (31),(32),(33) and (34) and Poincar´ e inequality we finally reach |bpod (un h,un h,en r)−bpod (un r,un r,en r)| ≤(1 +CC1,inf)∥en r∥2 0+C(C2 u,inf +C2 ld)∥∇ηn h∥2 0.(35) For the second term on the right-hand side of (29), we recall that bh(un h,un h,en r)=(un h·∇un h,en r)+1 2((∇ ·un h)un h,en r). bpod (un h,un h,en r)=(un h·∇un h,en r)+(en r·∇un h,un h)+((∇ ·un h)un h,en r). We notice that both bh(un h,un h,en r) and bpod (un h,un h,en r) share the term (un h· ∇un h,en r), and both also have the term ((∇ ·un h)un h,en r) but with a factor 1/2 in the second one, so that we have bh(un h,un h,en r)−bpod (un h,un h,en r)= −(en r·∇un h,un h)−1 2((∇ ·un h)un h,en r) Now, integrating by parts the first term on the right-hand side above we have bh(un h,un h,en r)−bpod (un h,un h,en r)=1 2(|un h|2,∇ ·en r)−1 2((∇ ·un h)un h,en r), and, since ∇ ·u=0, we can write bh(un h,un h,en r)−bpod (un h,un h,en r)=1 2(|un h|2,∇ ·en r)+1 2((∇ ·(un−un h))un h,en r),(36) For the second term on the right-hand side of (36) we have 1 2((∇ ·(un−un h))un h,en r)≤1 2∥un h∥∞∥∇ ·(un h−un)∥0∥en r∥0(37) ≤C2 u,inf µ 4∥∇ ·(un h−un)∥2 0+1 4µ∥en r∥2 0, 7 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 For the first term on the right-hand side of (36), we have 1 2(|un h|2,∇ ·en r)| ≤ 1 2(|un h|2−|un|2,∇ ·en r)|+ 1 2|(|un|2,∇ ·en r)| =1 2|(|un h|2−|un|2,∇ ·en r)|+ 1 2|((I−PQh)|un|2,∇ ·en r)|(38) ≤1 2∥un+un h∥∞∥un−un h∥0∥∇ ·en r∥0+1 2∥(I−PQh)|un|2∥0∥∇ ·en r∥0. Then, applying (16) 1 2|(|un h|2,∇ ·en r)| ≤ 1 2(∥un∥∞+Cu,inf)∥un−un h∥0∥∇ ·en r∥0 +1 2∥(I−PQh)|un|2∥0∥∇ ·en r∥0≤Cµ−1(∥un−un h∥2 0+∥(I−PQh)|un|2∥2 0) +µ 4∥∇ ·en r∥2 0,(39) where the generic constant Cabove depends on ∥u∥L∞(L∞)and Cu,inf. Inserting (35),(36),(37) and (39) into (29) we get |I| ≤ (1+CC1,inf +1 4µ)∥en r∥2 0+µ 4∥∇ ·en r∥2 0 +C(C2 u,inf +C2 ld)∥∇ηn h∥2 0+C2 u,inf µ 4∥∇ ·(un h−un)∥2 0 +Cµ−1(∥u−uh∥2 0+∥(I−PQh)|un|2∥2 0). Including into a generic constant Cthe dependence on Cu,inf and Cld of the third and fourth terms above we may write |I| ≤ (1+CC1,inf +1 4µ)∥en r∥2 0+µ 4∥∇ ·en r∥2 0 +C∥∇ηn h∥2 0+Cµ∥∇ ·(un h−un)∥2 0 +Cµ−1(∥u−uh∥2 0+∥(I−PQh)|un|2∥2 0).(40) We also have |I I| ≤ ν 2∥∇ηn h∥2 0+ν 2∥∇en r∥2 0,(41) and |III| ≤ ∥∇ ·ηn h∥2 0+µ 4∥∇ ·en r∥2 0.(42) Inserting (40),(41) and (42) into (28) and adding terms we get ∥en r∥2 0+ν n ∑ j=1 ∆t∥∇ej r∥2 0+µ n ∑ j=1 ∆t∥∇ ·ej r∥2 0 ≤ ∥e0 r∥2 0+ n ∑ j=1 ∆t(1+CC1,inf +1 4µ)∥ej r∥2 0+C(ν+µ) n ∑ j=1 ∆t∥∇ηj h∥2 0 +Cµ n ∑ j=1 ∆t∥∇ ·(uj h−uj)∥2 0+Cµ−1 n ∑ j=1 ∆t∥uj h−uj∥2 0 +Cµ−1 n ∑ j=1 ∆t∥(I−PQh)|uj|2∥2 0.(43) 8 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 Assuming ∆tCu:= ∆t(1 +CC1,inf +1 4µ)≤1 2(44) and applying Gronwall’s Lemma [15, Lemma 5.1] and (14),(15),(11) and (7) we obtain ∥en r∥2 0+ν n ∑ j=1 ∆t∥∇ej r∥2 0+µ n ∑ j=1 ∆t∥∇ ·ej r∥2 0 ≤e2T Cu(∥e0 r∥2 0+CT (ν+µ)∥Sv∥0 dv ∑ k=r+1 λk(45) +C(1 +µ−1T)C(u,p,l+1)2(h2l+(∆t)2)+Cµ−1h2lT∥|u|2∥L∞(Hl)). Remark 4.1. Let us observe that the error terms in the last line of (45) arise from using different discretization for the nonlinear term in the FOM and the POD methods. We also observe that these errors are of the same size as the error of the FOM method. More precisely, in view of (37) and (38), the extra errors have the size of µ∥∇·(un h−un)∥2 0, ∥u−uh∥2 0and ∥(I−PQh)|un|2∥2 0times several constants of the same size as the others appearing in the error analysis of the method. Although in general the dominant contribution to the extra error comes from µ∥∇ · (un h−un)∥2 0, thanks to the use of grad–div stabilization the error µ∥∇ ·(un h−un)∥2 0has the same rate of decay with has the L2 velocity error: ∥u−uh∥2 0. Theorem 4.2. Let ube the velocity in the Navier–Stokes equations (1), let urbe the grad–div POD stabilized approximation defined in (21), assume that the solution (u,p)of (1) is regular enough and condition (44) holds. Then, the error can be bounded as follows n ∑ j=1 ∆t∥uj r−uj∥2 0≤3T e2T Cu(∥e0 r∥2 0+CT (ν+µ)∥Sv∥0 dv ∑ k=r+1 λk +C(1 +µ−1T)C(u,p,l+1)2(h2l+(∆t)2) +Cµ−1h2lT∥|u|2∥L∞(Hl))+3T dv ∑ k=r+1 λk +3T C(u,p,l+1)2(h2l+(∆t)2).(46) Proof. Since ∑n j=1∆t∥ej r∥2 0≤Tmax1≤j≤n∥ej r∥2 0and n ∑ j=1 ∆t∥uj r−uj∥2 0≤3⎛ ⎝ n ∑ j=1 ∆t∥ej r∥2 0+ n ∑ j=1 ∆t∥Pv ruj h−uj h∥2 0 + n ∑ j=1 ∆t∥uj h−uj∥2 0⎞ ⎠, from (45),(14) and (11) we easily obtain (46).□ Remark 4.3. We have chosen, as in [1], to analyze the concrete case in which the divergence form is used for the FOM method and the EMAC form is applied for the POD method but the error analysis of any other combination of discretizations of the nonlinear terms with analogous properties could be carried out in a similar way. Remark 4.4. With the POD method defined in (21) we can obtain an approximation to the velocity but not to the pressure. In case an approximation to the pressure is also required one can use a supremizer pressure recovery 9 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 Fig. 5. Temporal evolution of drag coefficient computed with grad–div-ROM (21) (Skew, top; Emac, bottom) with µon =1 using r=8,12,16 velocity modes. Comparison with grad–div-FOM (8). Fig. 6. Temporal evolution of the strong (red) and weak (black) divergence of uh. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) 16 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 Fig. 7. Temporal evolution of drag coefficient computed with grad–div-ROM (21) with µon =0.1 using r=12 velocity modes. Comparison with grad–div-FOM (8). Fig. 8. Temporal evolution of drag coefficient computed with grad–div-ROM (21) (Skew, top; Emac, bottom) with µon =0.1 using r=8,12,16 velocity modes. Comparison with grad–div-FOM (8). 17 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 Fig. 9. Temporal evolution of kinetic energy (top) and drag coefficient (bottom) computed with ROM (21) (µ=0) using r=12 velocity modes. Comparison with FOM (8) (µ=0). equal to the value ∥∇ · uh∥0which is considerably big, explaining the big difference between the FOM values of the drag coefficient compared to the ROM values. To support the above explanation we now repeat Figs. 4 and 5with µon =0.1. In this case the difference |µon −µoff|is divided by around 10 and |µon −µoff|1/2is 0.3. In Figs. 7 and 8we have plotted the results. While in Figs. 4 and 5the error in the drag coefficients has the size of ∥∇ ·uh∥0, in Figs. 7 and 8this error has size 0.3 times ∥∇ ·uh∥0(the error is divided by 3.3 approximately respect to the previous experiment as can be observed). Since the error |µon −µoff|1/2∥∇ · uh∥0is smaller than before, in Fig. 8 we can now observe some improvement while increasing the number of modes. The improvement is bigger from r=8 to r=12 than from r=12 to r=16 since once the error reaches the quantity |µon −µoff|1/2∥∇ · uh∥0, no further reduction of the errors is expected. 5.3. Numerical results for ROM without grad–div and different discretization of the nonlinear term Finally, as in [1], we repeat the numerical experiments of Section 5.1 without grad–div stabilization, for which the error analysis performed in Section 4.2 holds. In particular, we use the FOM (8) with µ=0 to compute the snapshots for the construction of the POD modes. With POD velocity modes generated, the fully discrete ROM (21) with µ=0 is run with both skew-symmetric (as for the FOM) and EMAC (different from the FOM) formulation of the nonlinear term. Numerical results for kinetic energy and drag predictions using r=12 velocity modes are shown in Fig. 9, where we display a comparison of FOM (8) (µ=0) with the skew-symmetric form of the nonlinear term, and ROM (21) (µ=0) with both skew-symmetric and EMAC form of the nonlinear term. 18 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 Fig. 10. Temporal evolution of kinetic energy computed with ROM (21) (Skew, top; Emac, bottom; µ=0) using r=8,12,16 velocity modes. Comparison with FOM (8) (µ=0). From this figure, similarly to Section 5.1 and in agreement with [1], we see again that the ROM with the skew-symmetric form of the nonlinear term (as for the FOM) gives better results for both kinetic energy and drag predictions, in the sense that the good match that we observe using the same discretization for the nonlinear term as for the FOM (Skew) is slightly degraded using a different formulation of the nonlinear term (EMAC) in the online phase, as predicted by the theoretical study. Also in this context, we show numerical results for kinetic energy and drag predictions varying r, i.e. the number of velocity modes, in Figs. 10 and 11, respectively, where we display a comparison to FOM (8) (µ=0) with the skew-symmetric form of the nonlinear term. From these results, similarly to Section 5.1, we can again observe that better results are obtained when using the ROM with the skew-symmetric form of the nonlinear term (as for the FOM). In view of Remark 4.6, we can avoid the appearance of the term ∥∇ · (uh−u)∥0in the error analysis having instead the term ν−1/2∥uh−u∥0. Since ν−1/2≈31.6 in this experiment this extra term is still small and, consequently, we can observe a reduction of the ROM error when increasing rwith both nonlinear terms, although the reduction is larger when the same term is used in the FOM and ROM methods. We want to remark that, in case of using a coarser mesh to compute the FOM approximation and/or a bigger Reynolds number, the error term: min(ν−1/2∥uh−u∥0,∥∇·(uh−u)∥0) could dominate the behavior of the ROM approximation giving no improvement while increasing the number of modes, see the numerical experiments in [1]. 6. Conclusions As a conclusion, we can say that there are extra error terms in the case in which different discretizations are used for the nonlinear terms in the FOM and POD methods. Although these extra error terms lead in general to a less accurate method, in many examples, specially in the case in which grad–div stabilization is used, the extra error is small and the non consistent FOM–ROM method can still be used in practice. 19 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 Fig. 11. Temporal evolution of drag coefficient computed with ROM (21) (Skew, top; Emac, bottom; µ=0) using r=8,12,16 velocity modes. Comparison with FOM (8) (µ=0). In the case in which grad-did stabilization is added to both the FOM and POD methods error bounds with constants independent on inverse powers of the viscosity can be obtained. Comparing the case in which the same discretization of the nonlinear term is used for both FOM and POD methods we have bounded the added term in the error that has the size of the error of the FOM method (in the L2norm of the velocity and the L2norm of the divergence of the velocity). On the other hand, in the case in which no stabilization is added neither to the FOM nor to the POD method, we have carried out two different error analysis. In the first one, the error bounds depend on inverse powers of the viscosity and on the L2error of the velocity and the L2error of the divergence of the velocity of the FOM method plus the L2error of the projection of the square modulus of the velocity onto the FEM pressure space. In the second one, we are able to prove error bounds with constants independent on inverse powers of the viscosity but with the price of having bounds that depend on the H1norm of the FEM velocity error and the H1error of the projection of the square modulus of the velocity onto the FEM pressure space. Moreover, since the bounds for the plain Galerkin method depend on inverse powers of the viscosity the dependence on ν−1cannot completely be avoided also in this case. Overall, we have identified all the terms in the error bounds for different combinations of discretizations of the nonlinear term adding or not grad–div stabilization. Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. 20 B. Garc´ ıa-Archilla, J. Novo and S. Rubino Computer Methods in Applied Mechanics and Engineering 405 (2023) 115866 Data availability No data was used for the research described in the article. References [1] S. Ingimarson, L.G. Rebholz, T. Iliescu, Full and reduced order model consistency of the nonlinearity discretization in incompressible flows, Comput. Methods Appl. Mech. Engrg. 401 (part B) (2022) 16, Paper No. 115620. [2] J. Novo, S. Rubino, Error analysis of proper orthogonal decomposition stabilized methods for incompressible flows, SIAM J. Numer. Anal. 59 (1) (2021) 334–369. [3] S. Charnyi, T. Heister, M.A. Olshanskii, L.G. Rebholz, On conservation laws of Navier-Stokes Galerkin discretizations, J. Comput. Phys. 337 (2017) 289–308. [4] S. Charnyi, T. Heister, M.A. Olshanskii, L.G. Rebholz, Efficient discretizations for the EMAC formulation of the incompressible Navier-Stokes equations, Appl. Numer. Math. 141 (2019) 220–233. [5] R.A. Adams, Sobolev Spaces, in: Pure and Applied Mathematics, vol. 65, Academic Press, New York-London, 1975, [A subsidiary of Harcourt Brace Jovanovich, Publishers]. [6] V. John, Finite Element Methods for Incompressible Flow Problems, in: Springer Series in Computational Mathematics, vol. 51, Springer, Cham, 2016. [7] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, in: Classics in Applied Mathematics, vol. 40, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002, Reprint of the 1978 original [North-Holland, Amsterdam; MR0520174 (58 #25001)]. [8] F. Brezzi, R.S. Falk, Stability of higher-order Hood-Taylor methods, SIAM J. Numer. Anal. 28 (3) (1991) 581–590. [9] C. Taylor, P. Hood, A numerical solution of the Navier-Stokes equations using the finite element technique, Internat. J. Comput. Fluids 1 (1) (1973) 73–100. [10] J. de Frutos, B. García-Archilla, V. John, J. Novo, Analysis of the grad-div stabilization for the time-dependent Navier-Stokes equations with inf-sup stable finite elements, Adv. Comput. Math. 44 (1) (2018) 195–225. [11] K. Kunisch, S. Volkwein, Galerkin proper orthogonal decomposition methods for parabolic problems, Numer. Math. 90 (1) (2001) 117–148. [12] S. Ingimarson, An energy, momentum, and angular momentum conserving scheme for a regularization model of incompressible flow, J. Numer. Math. 30 (1) (2022) 1–22. [13] M.A. Olshanskii, L.G. Rebholz, Longer time accuracy for incompressible Navier-Stokes simulations with the EMAC formulation, Comput. Methods Appl. Mech. Engrg. 372 (2020) 113369, 17. [14] B. García-Archilla, V. John, J. Novo, On the convergence order of the finite element error in the kinetic energy for high Reynolds number incompressible flows, Comput. Methods Appl. Mech. Engrg. 385 (2021) 54, Paper No. 114032. [15] J.G. Heywood, R. Rannacher, Finite-element approximation of the nonstationary Navier-Stokes problem. IV. Error analysis for second-order time discretization, SIAM J. Numer. Anal. 27 (2) (1990) 353–384. [16] K. Kean, M. Schneier, Error analysis of supremizer pressure recovery for POD based reduced-order models of the time-dependent Navier-Stokes equations, SIAM J. Numer. Anal. 58 (4) (2020) 2235–2264. [17] L.R. Scott, M. Vogelius, Conforming finite element methods for incompressible and nearly incompressible continua, in: Large-Scale Computations in Fluid Mechanics, Part 2 (La Jolla, Calif., 1983), in: Lectures in Appl. Math., vol. 22, Amer. Math. Soc., Providence, RI, 1985, pp. 221–244. [18] R.S. Falk, M. Neilan, Stokes complexes and the construction of stable finite elements with pointwise mass conservation, SIAM J. Numer. Anal. 51 (2) (2013) 1308–1326. [19] M.A. Case, V.J. Ervin, A. Linke, L.G. Rebholz, A connection between Scott-Vogelius and grad-div stabilized Taylor-Hood FE approximations of the Navier-Stokes equations, SIAM J. Numer. Anal. 49 (4) (2011) 1461–1481. [20] M. Schäfer, S. Turek, Benchmark computations of laminar flow around a cylinder, in: E.H. Hirschel (Ed.), Flow Simulation with High-Performance Computers II, in: Notes on Numerical Fluid Mechanics, vol. 48, Vieweg, 1996, pp. 547–566. [21] F. Hecht, New development in freefem++, J. Numer. Math. 20 (3–4) (2012) 251–265. [22] N. Ahmed, S. Rubino, Numerical comparisons of finite element stabilized methods for a 2D vortex dynamics simulation at high Reynolds number, Comput. Methods Appl. Mech. Engrg. 349 (2019) 191–212. [23] S. Giere, T. Iliescu, V. John, D. Wells, SUPG reduced order models for convection-dominated convection–diffusion-reaction equations, Comput. Methods Appl. Mech. Engrg. 289 (2015) 454–474. [24] V. John, Reference values for drag and lift of a two-dimensional time-dependent flow around a cylinder, Internat. J. Numer. Methods Fluids 44 (2004) 777–788. [25] B. García-Archilla, J. Novo, S. Rubino, Error analysis of proper orthogonal decomposition data assimilation schemes with grad-div stabilization for the Navier-Stokes equations, J. Comput. Appl. Math. 411 (2022) 30, Paper No. 114246. [26] C. Zerfas, L.G. Rebholz, M. Schneier, T. Iliescu, Continuous data assimilation reduced order models of fluid flow, Comput. Methods Appl. Mech. Engrg. 357 (2019) 112596, 18. 21