Desarrollo y análisis de un modelo de orden reducido POD con viscosidad artificial. Aplicaciones en mecánica de fluidos

Trabajo Fin de M´aster: Desarrollo y an´alisis de un modelo de orden reducido POD con viscosidad artificial. Aplicaciones en mec´anica de fluidos Carlos N´u˜nez Fern´andez M´aster Univesitario en Matem´aticas Tutorizado por: Samuele Rubino Febrero 2022 DPTO. ECUACIONES DIFERENCIALES Y AN´ ALISIS NUM´ ERICO - US 2 Resumen La simulaci´on de flujos de flu´ıdos con altos n´umeros de Reynolds (r´egimen de convecci´on dominante, turbulento) sigue siendo un ´area de investigaci´on muy activa. La resoluci´on precisa y eficiente al mismo tiempo de las ecuaciones en derivadas parciales subyacentes puede suponer un verdadero desaf´ıo. Por esta raz´on, la finalidad de este trabajo es proponer Modelos de Orden Reducido (MOR), de tipo Proper Orthogonal Descomposition (POD), para tratar este problema. En particular, en este trabajo se expone un estudio te´orico del m´etodo POD y se propone un modelo de orden reducido POD con viscosidad artifical. Este m´etodo se aplicar´a en primera instancia a la ecuaci´on de Burgers con coeficientes de difusi´on muy peque˜nos. En este contexto, para el m´etodo propuesto (AV-POD-G-ROM) se derivan las estimaciones de error correspondientes y se realizan simulaciones num´ericas. En segundo lugar, se extiende este m´etodo a las ecuaciones de Navier- Stokes con viscosidades peque˜nas, dando lugar a un modelo de orden reducido POD de tipo Smagorinsky (S-POD-G-ROM). En este contexto, se realizan simulaciones num´ericas del flujo no estacionario alrededor de un obst´aculo. 3 4 Abstract The numerical simulation of fluid flows at high Reynolds numbers (dominant convection, turbulent regime) is a very active research area. The precise and efficient resolution of the corresponding partial derivative equations can be a real challenge. For this reason, the goal of this work is to propose Reduced Order Models (ROMs), based on Proper Orthogonal Decomposition (POD), to treat this problem. In particular, in this work we show a theoretical study of the POD method and propose a POD reduced order model with artificial viscosity. This method will be applied in first instance to the Burgers equation with small diffusion coefficients. In this context, for the proposed method (AV-POD-G-ROM) the corresponding error estimates are derived and numerical simulations will be carried out. Secondly, this method will extend to the Navier-Stokes equations with small viscosity, leading to a Smagorinsky POD reduced order model (S-POD-G-ROM). In this context, numerical simulations of unsteady flow around an obstacle will be carried out. 5 6 Contents Resumen 3 Abstract 5 1 Introduction 13 2 Notation and preliminary results 15 2.1 Notation................................. 15 2.2 Preliminaryresults ........................... 16 3 The POD method in Rmspace 19 3.1 Singular Value Decomposition . . . . . . . . . . . . . . . . . . . . . 19 3.2 The relation between POD and SVD . . . . . . . . . . . . . . . . . 20 3.3 Application to Image Compression . . . . . . . . . . . . . . . . . . 24 3.4 The POD method with weighted inner product . . . . . . . . . . . 26 3.5 Treatment of non-linear and non-affine terms within the POD method 28 3.5.1 Discrete Empirical Interpolation Method . . . . . . . . . . . 28 3.5.2 Radial Basis Functions . . . . . . . . . . . . . . . . . . . . . 30 4 The POD for the Burgers equation 33 4.1 AV-POD-G-ROM............................ 34 4.2 Errorestimates ............................. 36 4.3 Practical Implementation and Numerical Results . . . . . . . . . . . 41 5 The POD for the Navier-Stokes equations 49 5.1 S-POD-G-ROM............................. 52 5.2 NumericalStudies............................ 53 5.2.1 Setup for numerical simulations . . . . . . . . . . . . . . . . 53 5.2.2 FOM and POD modes . . . . . . . . . . . . . . . . . . . . . 55 5.2.3 Adaptive in time algorithm for the Smagorinsky coefficient Cs................................ 58 5.2.4 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . 58 7 8CONTENTS 6 Conclusion and outlook 63 6.1 Conclusion................................ 63 6.2 Outlook ................................. 63 6.2.1 NSE error estimates . . . . . . . . . . . . . . . . . . . . . . 63 6.2.2 Pressure recovery and turbulent flows . . . . . . . . . . . . . 64 A Ongoing research 65 Bibliography 68 List of Figures 3.1 Approximation of the original images with lower rank . . . . . . . . 25 3.2 Singular values of red, blue and green. . . . . . . . . . . . . . . . . 25 3.3 Energy ratio of red, blue and green. . . . . . . . . . . . . . . . . . 26 4.1 PODmodes............................... 42 4.2 DNSsolution............................... 42 4.3 POD-G-ROM solution. . . . . . . . . . . . . . . . . . . . . . . . . . 43 4.4 AV-POD-G-ROM solution in space-time with ˜r=20 ........ 43 4.5 AV-POD-G-ROM solution in space-time with ˜r=40 ........ 44 4.6 AV-POD-G-ROM solution in space-time with ˜r=60 ........ 44 4.7 Solution curves of DNS, POD-G-ROM and AV-POD-G-ROM (˜r= 20) at final time T= 1.......................... 45 4.8 Solution curves of DNS, POD-G-ROM and AV-POD-G-ROM (˜r= 40) at final time T= 1.......................... 46 4.9 Solution curves of DNS, POD-G-ROM and AV-POD-G-ROM (˜r= 60) at final time T= 1.......................... 46 5.1 Computationalgrid........................... 54 5.2 Final velocity FOM. . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.3 Temporal evolution of drag and lift coefficients. . . . . . . . . . . . 56 5.4 Temporal evolution of kinetic energy. . . . . . . . . . . . . . . . . 57 5.5 Decay of the eigenvalues and captured energy. . . . . . . . . . . . . 57 5.6 First four POD modes of velocity. . . . . . . . . . . . . . . . . . . . 59 5.7 Time evolution of the kinetic energy for FOM, POD-G-ROM, S- POD-G-ROM and Adapt-S-POD-G-ROM for r=7. ........ 60 5.8 Temporal evolution of absolute error in kinetic energy of POD-G- ROM, S-POD-G-ROM and Adapt-S-POD-G-ROM with respect to FOM on the left and the value of Cscoefficient in the adaptive method on the right using r=7. ................... 60 9 16 CHAPTER 2. NOTATION AND PRELIMINARY RESULTS 4. Divergence of a vector function v:Rd→Rd: ∇·v= d X k=1 ∂kvk. 5. Divergence of a tensor-valued function σ:Rd→Rd×d: ∇·σ=  Pd k=1 ∂kσ1k . . . Pd k=1 ∂kσdk   . 6. Laplace operator of a scalar-valued function v:Rd→R: ∆v=∇·∇v= d X k=1 ∂kkv. 7. Convection term of a vector-valued function v:Rd→Rd, velocity of the Navier-Stokes equation: (v·∇)v=   v·∇v1 . . . v·∇vd   =  Pd k=1(∂kv1)vk . . . Pd k=1(∂kvd)vk   . 2.2 Preliminary results In this dissertation, we work with several function spaces and we use its associated norms, so we are going to introduce them. Let 1 ≤p≤ ∞,the Lebesgue spaces are defined as: Lp(Ω) = {f: Ω →R/f measurable, kfkLp(Ω) <∞}. The norm for a function u: Ω →Rof Lp(Ω) is kukLp(Ω) = (ZΩ|u(x)|pdx)1/p with 1 ≤p < ∞, kukL∞(Ω) = ess sup x∈Ω|u(x)|. Lebesgue spaces are Banach spaces and for p= 2 we get the Hilbert space L2(Ω) with inner product: 2.2. PRELIMINARY RESULTS 17 (u, v)L2(Ω) =ZΩ u(x)v(x)dx. Normally, we will use the following notation for the L2−inner product (·,·) := (·,·)L2(Ω). Sobolev spaces too play an important role in the analysis of PDEs. In this work, we use: H1(Ω) = {f∈L2(Ω)/∇f∈[L2(Ω)]d}, and when working with homogeneous Dirichlet boundary conditions, we consider: H1 0(Ω) = {f∈H1(Ω)/f|∂Ω= 0}. This is a closed linear subspace of H1and thus a Hilbert space endowed with the H1−norm (kukH1).Thanks to Poincar`e inequality, the H1−norm is equivalent on H1 0to the norm kukH1 0=k∇ukL2. Finally, we introduce the Bochner spaces, needed for the analysis of timedependent PDEs [9]. For 1 ≤p≤ ∞,the Bochner space is defined as: Lp((0, T),Ω) = {f: (0, T)×Ω→R/f Bochner measurable,kfkLp((0,T ),Ω) <∞}, where for a function u: (0, T)×Ω→R,its associated norm is: kukLp((0,T),Ω) = (ZT 0ku(t)kp Lp(Ω)dt)1/p with 1 ≤p < ∞, kukL∞((0,T),Ω) = ess sup t∈(0,T )ku(t)kL∞(Ω). 18 CHAPTER 2. NOTATION AND PRELIMINARY RESULTS Chapter 3 The POD method in Rmspace In this chapter, we introduce the Proper Orthogonal Decomposition (POD) method in Rm. This method seeks a proper orthonormal basis called POD basis {ϕi}r i=1 of rank r, which preserve the essential information of nvectors y1, . . . , yn∈Rm called snapshots, with r≤min{m, n}.Snapshots constitute a high dimensional data set that normally comes from a Direct Numerical Simulations(DNS) [10, 27]. The POD method is formulated as constrained optimization problem: max ϕ1,...,ϕr∈Rm r X i=1 n X j=1 (yj, ϕi)Rmsubject to (ϕi, ϕj)Rm=δij 1≤i, j ≤r, (3.1) where δij is the Kronecker delta and (a, b)Rm=Pm i=1 aibi. In this chapter, we are going to introduce the Singular Value Descomposition (SVD) and the connection with the POD, and some properties of the POD basis. This part follows [30]. Afterward, we consider a widely used application of the SVD method, such as image compression. We introduce the POD method with weighted inner product, useful when applying it to PDEs. Finally, we consider the treatment of non-linear and non-affine terms (with respect to the parameter) within the POD. 3.1 Singular Value Decomposition Let S= (y1|···|yn)∈Rm×nbe a real matrix of rank d≤min(m, n),usually called the snapshots matrix. The Singular Value Decomposition (SVD) of Sguarantees that there exists orthogonal matrices: Φ = [ϕ1|···|ϕm]∈Rm×mand V= [v1|. . . |vn]∈Rn×n(3.2) 19 20 CHAPTER 3. THE POD METHOD IN RMSPACE such that ΦTSV = Σ := D0 0 0 ∈Rm×n,(3.3) where D= diag(σ1, . . . , σd)∈Rd×d,with σ1> . . . > σd>0.(see [11]). Moreover the vectors {ϕi}d i=1 and {vi}d i=1 satisfy Sϕi=σiviand STvi=σiϕifor i= 1, . . . , d, (3.4) so STSϕi=σiSTvi=σ2 iϕi∀i= 1, . . . , d, (3.5) SSTϕi=σiSϕi=σ2 ivi∀i= 1, . . . , d. (3.6) Therefore, we obtain that {ϕi}d i=1 and {vi}d i=1 are the eigenvectors for SSTand STSand their eigenvalues are λi=σ2 i>0i= 1, . . . , d. From (3.3) and using that Φ and Vare orthogonal, we obtain that S= ΦΣVT.(3.7) Simplifying (3.7) we note S= ΦdD(Vd)T= ΦdBd,(3.8) where Φd∈Rd×n, V d∈Rd×nand Bd=D(Vd)T∈Rd×n.The expression (3.8) says that the column space of S can be represented in terms of the dlinearly independent columns of Φdand we have that sj=Pd i=1 Bd ijϕd i=Pd i=1(D(Vd)T)ijϕi=Pd i=1((Φd)TΦdD(Vd)T)ijϕj= (3.8) =Pd i=1((Φd)TS)ijϕj=Pd i=1(Pm k=1 ϕd kiskj)ϕj=Pd i=1(ϕi, sj)Rmϕi. Therefore, the column vector yjof Scan be expressed as linear combination of {ϕi}d i=1 as: yj= d X i=1 (ϕi, yj)Rmui∀j= 1, . . . , n. (3.9) 3.2 The relation between POD and SVD Now, we are going to connect the POD with the SVD. Firstly, we want to compute iteratively an orthonormal basis, which approximates the data set {yj}n j=1. From this computation, we obtain the vectors {ϕi}r i=1 that will be called the POD basis of rank r. The first POD mode solves the following maximization problem max ˜ϕ1∈Rm n X j=1 |(yj,˜ϕ1)Rm|2s.t. k˜ϕ1k2 Rm= 1.(3.10) 3.2. THE RELATION BETWEEN POD AND SVD 21 Thanks to Lagrange formalism, the first lest singular vector ϕ1solves (3.10). Following this procedure, the second POD mode is obtained by solving the following maximization problem max ˜ϕ2∈Rm n X j=1 |(yj,˜ϕ2)Rm|2s.t. k˜ϕ2k2 Rm= 1,(ϕ1,˜ϕ2)Rm= 0,(3.11) and the solution is given by the left singular vector ϕ2.According this inductive procedure, we get the following theorem: Theorem 1. Let S= [y1, . . . , yn]∈Rm×nwith rank d≤min(m, n).Moreover let S= ΦΣVTbe the SVD described in (3.7). Then for 1≤r≤dthe optimization problem max ˜ϕ1,..., ˜ϕr∈Rm r X i=1 n X j=1 |(yj,˜ϕi)Rm|2s.t. ( ˜ϕi,˜ϕj)Rm=δij∀1≤i, j ≤r, (3.12) is being solved by the left singular vectors {ϕi}r i=1 and it holds that arg max (3.12) = r X i=1 σ2 i= r X i=1 λi.(3.13) Proof. The proof of Theorem (1) can be found in [[30],theorem 1.1]. Having made the connection between the POD basis and the SVD through Theorem (1), we are going to show that the POD basis is optimal in the following sense: Theorem 2. (Optimality of the POD basis) Let all hypotheses of Theorem (1) be satisfied. Suppose that b Φd∈Rm×ddenotes a matrix with pairwise orthonormal vectors bϕiand that the expansion of the columns of S in the basis {bϕi}d i=1 is given by S=b ΦdCd, where Cd ij = ( bϕi, yj)Rmfor 1≤i≤d, 1≤j≤n. (3.14) Then for every r∈ {1, . . . , d}we have kS−ΦrBrkF≤ kS−b ΦrCrkF.(3.15) In (3.15),k·kFdenotes the Frobenius norm given by kAkF=v u u t m X i=1 n X j=1 |A2 ij|=ptrace(AtA)for A∈Rm×n,(3.16) the matrix Φrdenotes the first r≤dcolumns of Φ, Brthe first r rows of B and similarly for b Φrand Cr. 22 CHAPTER 3. THE POD METHOD IN RMSPACE To prove Theorem (2), we need the following result: Lemma 1. Let A∈Rm×nand Qan orthogonal matrix of order m. Then: kQAkF=kAkF.(3.17) Proof. kQAk2 F= trace((QA)T(QA)) = trace(ATQTQA) = trace(ATA) = kAkF.(3.18) Now we are going to prove Theorem (2). Proof. Using the Lemma (1) kY−b ΦrCrk2 F=kb Φd(Cd−Cr 0)k2 F=kCd−Cr 0k2 F= d X i=r+1 n X j=1 |Cr ij|2,(3.19) where Cr 0∈Rr×nis obtained by replacing the last d−rrows of C∈Rr×nby 0. Similarly we have kY−ΦrBrk2 F=kΦk(Bd−Br 0)k2 F=kBd−Br 0k2 F=Pd i=r+1 Pn j=1 |Bd ij|2 =Pd i=r+1 Pn j=1 |(yj, ϕi)|=Pd i=r+1 Pn j=1((yj, ϕi)yj, ϕi) = Pd i=r+1(SSTϕi, ϕi) = Pd i=r+1 σ2 i. (3.20) We know that, by Theorem (1), ϕ1, . . . , ϕdsolve problem (3.12). Now, using (3.20), we obtain: kSk2 F=kb ΦdCdk2 F=kCdk2 F= d X i=1 n X j=1 |Cr ij|2,(3.21) and kSk2 F=kΦdBdk2 F=kBdk2 F= d X i=1 n X j=1 |Br ij|2= r X i=1 σ2.(3.22) Therefore kS−ΦrBrk2 F=Pd i=r+1 σ2 i=Pd i=1 σ2 i−Pr i=1 σ2 i=kSk2 F−Pr i=1 Pn j=1 |(yj, ϕi)|2 ≤ kSk2 F−Pr i=1 Pn j=1 |(yj,bϕi)|2=Pd i=1 Pn j=1 |Cd ij|2−Pr i=1 Pn j=1 |Cd ij|2 =Pd i=r+1 Pn j=1 |Cd ij|2=kS−b ΦrCrk2 F. (3.23) 3.2. THE RELATION BETWEEN POD AND SVD 23 We note that kS−b ΦrCrk2 F= n X j=1 m X i=1 |Sij − r X k=1 (cϕk, yj)b Φr ik|2= n X j=1 kyj− r X k=1 (yj,cϕk)cϕkk2,(3.24) and kS−ΦrBrk2 F= n X j=1 kyj− r X k=1 (yj, ϕk)ϕkk2.(3.25) Hence, by Theorem (2) n X j=1 kyj− r X k=1 (yj, ϕk)ϕkk2≤ n X j=1 kyj− r X k=1 (yj,cϕk)cϕkk2.(3.26) Then, we can obtain the POD basis of rank rusing the Theorem (2) by solving the following optimization problem: min ˜ϕ1,..., ˜ϕr∈Rm n X j=1 kyj− r X i=1 (yj,˜ϕi) ˜ϕik2s.t. ( ˜ϕi,˜ϕj) = δij for 1 ≤i, j ≤r. (3.27) In practice, computing the SVD of S is computationally expensive, since S∈ Rm×nis a dense matrix and mor nare very large. However, we can transform the search of the POD basis into an eigenvalues problem of size m×mor n×n, which is useful if m << n or m >> n. Hence, we are going to distinguish two cases: 1. If n < m : We compute the eigenvectors φ1, . . . , φr∈Rnby solving the symmetric n×n eigenvalues problem: STSφi=λiφifor i= 1, . . . , r. (3.28) and using (3.4), ϕi=1 √λi Sφifor i= 1, . . . , r. (3.29) 2. If n > m: We can obtain the POD basis by solving the m×meigenvalues problem: SSTϕi=λiϕifor i= 1, . . . , r. (3.30) This method is usually called the method of snapshots for historical reasons [26]. 24 CHAPTER 3. THE POD METHOD IN RMSPACE Finally, we mention that for the choice of rin practical problems one could use the following formula: (r) := Pr i=1 λi Pd i=1 λi≥δ,(3.31) where (d) is called the energy ratio and δ∈(0,1] is a fixed threshold. Since we want to capture most of the energy system by the POD modes, usually we choose δclose to 1. 3.3 Application to Image Compression In this section, we show how helpful is the SVD in real life, showing one of its most used application that is image compression. Image compression [7] is based on: given an image S∈Rm×nof rank d, we are looking for a matrix ˜ S∈Rm×nwith rank r << d such that kS−˜ Skis minimal. Then by the Eckart-Young theorem [12], we have that: min rank(˜ S)=rkS−˜ Sk=kS−S(r)k=(σr+1 for the 2-norm. qPd i=r+1 σ2 ifor the Frobenius-norm. (3.32) The results (3.32) means that truncated matrix from the SVD: S(r):= Φ Σ(r)0 0 0 VT= r X i=1 σiϕivT i,(3.33) with Σ(r):= diag(Σ11,...,Σrr) solves the minimization problem. Consequently, we only need to compute a truncated SVD with a good rank (high content of information) in order to compress an image. We use a colour image of the Faculty of Mathematics (University of Seville), which is 600 pixels wide and 900 pixels high. Each pixel is divided in three integer values corresponding to the intensity of red, blue and green of the colour. In Figure (3.1), we can see the original image and its SVD approximation of rank 5, 150 and 500. 3.3. APPLICATION TO IMAGE COMPRESSION 25 (a) Original image (b) SVD Approximation of rank 5 (c) SVD Approximation of rank 150 (d) SVD Approximation of rank 500 Figure 3.1: Approximation of the original images with lower rank The matrix representation of the image has full rank, i.e., it has rank 600. First of all, we have to divide the matrix of the original image into three matrices corresponding to each colour. Then, we have to perform an SVD approximation for each colour, therefore the decayment of its singular values are shown in Figure (3.2). Figure 3.2: Singular values of red, blue and green. From Figures (3.2), we observe a fast decay of the corresponding singular values, so that we expect a good approximation for relatively small values of r. 32 CHAPTER 3. THE POD METHOD IN RMSPACE (c) This corresponds to the following linear system: Aξ LwL=YL,(3.44) where (Aξ L)ij =ξL,j(kxi µ,t −xj µ,tkL2).(3.45) The weights wLcan be thus computed in the offline phase, and then stored to be used in the online phase. (d) In the online stage, we have z∗and the goal is to compute α(z∗) = [αi(z∗)]˜r i=1 which it is done by α(z∗)≈Gi(z∗) = ˜ M X j=1 wi,jξi,j(kz∗−xj µ,tkL2(Rq+1))i= 1,2,...,˜r. (3.46) This technique has been used in Appendix (A) to approximate the eddy viscosity. Chapter 4 The POD for the Burgers equation The Burgers equations with initial and boundary conditions is given by:    ut−νuxx +uux=fin Ω ×(0, T), u(x, 0) = u0in Ω, u(x, t) = gon ∂Ω×(0, T). (4.1) where νis the diffusion parameter, fis the forcing term, Ω is the 1D computational domain and Tis the final time. We are going to consider the function space X=H1 0(Ω) and assume that, without loss of generality, g≡0 and f≡0.Then, the weak formulation of the Burgers equation (4.1) reads: (ut, v)L2+ν(ux, vx)L2+ (uux, v)L2= 0 ∀v∈X. (4.2) To derive the standard POD-Galerkin ROM applied to the Burgers equation, one has to solve first the fully discretized version of (4.2) for a certain set of time instances. In this way, we collect the snapshots, and we perform the POD method over the set of snapshots to compute the POD basis in L2(see Section 3.4). The Galerkin projection-based POD-ROM uses both Galerkin truncation and Galerkin projection. The former yields an approximation of the solution by a linear combination of the truncated POD basis: u(x, t)≈ur(x, t) = r X j=1 aj(t)ϕj(x)∈Xr,(4.3) where Xr= span{ϕ1, . . . , ϕr}is the space for the POD setting and {ai(t)}r i=1 are the time coefficients. Replacing the solution uwith urin the Burgers equations (4.2), using the Galerkin method, and projecting the resulted equations onto the space Xr, one 33 34 CHAPTER 4. THE POD FOR THE BURGERS EQUATION obtains the standard POD-Galerkin ROM (POD-G-ROM) for the Burgers equations. To obtain its time discretization, we apply the backward Euler method and denoting ¯ ∂uk r=uk r−uk−1 r ∆t, the problem reads: (¯ ∂uk r, ϕ)L2+ν((uk r)x,(ϕ)x)L2+ (uk r(uk r)x, ϕ)L2= 0 ∀ϕ∈Xr.(4.4) The POD-G-ROM (4.4) leads to the following autonomous system for the coefficients a(t) : · a(t) = Ba +aTCa, (4.5) where Band Ccorrespond to the linear and quadratic terms in the numerical discretization of the Burgers equation (4.4), respectively. The initial conditions are obtained by the projection: · ai(0) = (u(·,0), ϕi)L2i= 1, . . . , r. The finite dimensional system (4.5) can be written componentwise as follow: · ai(t) = r X m=1 Bimam(t) + r X n=1 r X m=1 Cimnam(t)an(t). where Bim =−((ϕm)x,(ϕi)x)L2, Cimn =−(ϕm(ϕn)x, ϕi)L2. 4.1 AV-POD-G-ROM In this section we are going to introduce an artificial viscosity POD-ROM (AVPOD-G-ROM). The reason to introduce this model is that for convection-dominated non-linear problems, the standard POD-G-ROM is not accurate. Thus, we introduce a new term of artificial viscosity (AV), which try to replicate in this framework the effect of the Smagorinsky turbulence model for the Navier-Stokes equations [32]. We emphasize that, although the Burgers equations is considered as simplification of NSE and small diffusion parameters νare corresponding to high Reynolds numbers in NSE, in this context we are not representing the real turbulence. Therefore, the role of closure methods changes correspondingly. Instead of improving physical accuracy (dissipating energy according to the well known concept of energy cascades), they increase the numerical stability (in testing problems 4.1. AV-POD-G-ROM 35 showing shock-like behavior, for instance). Of course, once we fully understand the behavior of this kind of POD closure model in this simplified setting, we will analyze in the NSE setting (see (6)). Therefore, the new weak formulation of (4.2) is: (ut, v)L2+ν(ux, vx)L2+ (uux, v)L2+c(|ux|ux, vx)=0 ∀v∈X, (4.6) where cis a positive coefficient. Now, we repeat the same process as in the previous section, so we replace (4.3) into (4.6) and obtain the AV-POD-G-ROM: (¯ ∂uk r, ϕ)L2+ν((uk r)x,(ϕ)x)L2+ (uk r(uk r)x, ϕ)L2+c(|uk r|uk r, φx)L2= 0 ∀ϕ∈Xr. (4.7) In order to obtain an autonomous system with respect to the time coefficients a(t),we have to treat the non-linear term (|(ud)x|(ud)x, ϕx).In this case we are going to do an approximation by the DEIM (see Section (3.5.1)). Up to our knowledge, it is the first time that the DEIM is applied in such kind of model to reduce its computational complexity, as we will show in Section (4.3). Then, if we apply the DEIM method on |ux|we obtain: |ux|= ˜r X j=1 αj(t) ˜ϕj, where ˜ X˜r=span{˜ϕ1,..., ˜ϕ˜r}and αj(t) are the time coefficients. Consequently we obtain the following autonomous system for the coefficients a(t) : · a(t) = Ba +aTCa +αTDa. (4.8) where B, C and Dcorrespond to the linear and the two quadratic terms in the numerical discretization of the AV Burgers equation (4.7), respectively. The initial conditions are obtained by the projection: · ai(0) = (u(·,0), ϕi)L2i= 1, . . . , r. The finite dimensional system (4.8) can be written componentwise as follow: · ai(t) = r X m=1 Bimam(t) + r X n=1 r X m=1 Cimnam(t)an(t) + ˜r X n=1 r X m=1 Dimnam(t)αn(t), where 36 CHAPTER 4. THE POD FOR THE BURGERS EQUATION Bim =−((ϕm)x,(ϕi)x)L2, Cimn =−(ϕm(ϕn)x, ϕi)L2, Dimn =−(c( ˜ϕm)x(ϕn)x,(ϕi)x)L2. 4.2 Error estimates In this section, our aim is to prove estimates for the time discretization error 1 MPM k=1 kuk r−u(tk)k2 L2,where uk ris the solution of (4.7) and u(tk) is the solution of (4.6). For it, we follow [19] and [32]. However, in this section we only focus on time discretization and POD truncation, for finite element spatial discretization the reader is referred to [20] and [15]. First of all, we are going to introduce some notations and results that we need to prove the main result. For clarity and simplicity, let L2, H1 0, W1,3, W1,3 0denote the usual Sobolev spaces on Ω and (·,·) = (·,·)L2. For all ϕ, ψ ∈H1 0,we define a(ϕ, ψ) := ν(ϕx, ψx)∀ϕ, ψ ∈H1 0,(4.9) (F(ϕ), ψ) := (ϕϕx, ψ).(4.10) For ϕ, ψ ∈W1,3and c > 0, (G(ϕ), ψ) := c(|ϕx|ϕx, ψx).(4.11) We will denote by Ck, k ∈Nthe different constants that will appear independent of the number of the POD modes (r) and the number of snapshots (M). Now, we show several results, that will be used throughout this section. The first one: Lemma 2 (Poincar`e inequality).∃α > 0constant such that kϕkL2≤α|ϕ|H1∀ϕ∈H1 0,(4.12) where |·|H1is the H1−seminorm. Lemma 3. The bilinear form a(·,·)is continuous and coercive, this means that ∃β, κ such that |a(ϕ, ψ)| ≤ βkϕkH1kψkH1∀ϕ, ψ ∈H1 0,(4.13) |a(ϕ, ϕ)| ≤ κkϕk2 H1∀ϕ∈H1 0.(4.14) 4.2. ERROR ESTIMATES 37 Lemma 4 (Strong monotonicity and Lipschitz continuity of G).For all ϕ1, ϕ2, ψ ∈ W1,3,∃Cconstant, such that (G(ϕ1)−G(ϕ2), ϕ1−ϕ2)≥Ck(ϕ1−ϕ2)xk3 L3,(4.15) (G(ϕ1)−G(ϕ2), ψ)≤C˜ Mk(ϕ1−ϕ2)xkL3kψxkL3,(4.16) where ˜ M:= max{k(ϕ1)xkL3,k(ϕ2)xkL3}. Definition 1. The projection PR:H1 0→Xris defined as follows: for all ψ∈Xr a(PRu, ψ) = a(u, ψ)∀u∈H1 0.(4.17) Lemma 5. If ∆tis small enough, we have the following approximation for PR([19]): 1 M M X k=1 ku(k)−PRu(tk)k2 H1≤C1 d X j=r+1 λj,(4.18) 1 M M X k=1 k¯ ∂u(tk)−PR∂¯u(tk))k2 H1≤C2 d X j=r+1 λj.(4.19) We also need the following results: Lemma 6. ∃C3, C4such that for all v∈Xit is verified kvkL3≤C3kvkL2kvkL2≤C4kvkL3.(4.20) Proof. We note that {ϕ1, . . . , ϕd}is an orthonormal basis of X, so if v∈X, then ∃viscalars such that vi= (v, ϕi) and v=Pd i=1 viϕi.Consequently kvkL3=k d X i=1 viϕikL3≤ d X i=1 |vi|kϕikL3≤( max i=1,...,r kϕikL3) d X i=1 |vi|.(4.21) Now we need the following auxiliar Lemma (7) whose proof is in ([5]): Lemma 7 (auxiliar).Let {x1, . . . , xn}be a linearly independent set of vectors in a normed space X (of any dimension). Then there is a number C > 0such that for every choice of scalars α1, . . . , αnwe have kα1x1+. . . +αnxnk ≥ C(|α1|+···+|αn|).(4.22) Applying the auxiliar lemma (7) into (4.21) we obtain the first inequality. The other one has the same proof as the first. 38 CHAPTER 4. THE POD FOR THE BURGERS EQUATION Lemma 8. Assume that the solution of (4.6) satisfies ux∈H1 0.There exists constants C5, C6and C7, such that, for all k= 1, . . . , m, the following inequalities hold kux(tk)kL2≤C5,kux(tk)kL3≤C6,kPRux(tk)kL2≤C7.(4.23) Proof. The first two inequalities are easily proved by letting ϕ=uin (4.6), using (F(u), u) = 0,Cauchy-Schwarz inequality, the strong monotonicity of G (4.15) and adjusting constants. Therefore, we want to prove the last inequality, so first of all, we remind the definition of the projection (4.17): a(PRux, ψ) = a(ux, ψ)∀ψ∈Xr.(4.24) If we let ψ=PRux,we obtain: a(PRux, PRux) = a(ux, PRux).(4.25) Next, we are going to bound separately both terms in (4.25). Firstly, using the continuity of a(·,·) : a(PRux, ux)≤βkPRuxkL2kuxkL2.(4.26) Secondly, using the coercivity of a(·,·) : a(PRux, PRux)≥κkPRuxk2 L2.(4.27) Then, if we join (4.27),(4.26) and (4.25): κkPRuxk2 L2≤a(PRux, PRux) = a(ux, PRux)≤βkPRuxkL2kuxkL2,(4.28) kPRuxkL2≤β κkuxkL2≤C7.(4.29) We are now ready to prove the main theorem: Theorem 4. Let u be the solution of (4.6) and {uk r}M k=1 be the solution of (4.7). Assume that ux∈H1 0and utt ∈L2(0, T;L2).If ∆tis small enough, then there exists a constant C > 0, such that 1 M M X k=1 kuk r−u(tk)kL2≤C(ku0−PRu0k2+ d X j=r+1 λj+ ∆t2).(4.30) 4.2. ERROR ESTIMATES 39 Proof. The proof follows the same ideas and procedure as the Theorem 7 and 11 in ([19]) and Theorem 1 in section 4 of [32]. For that reason, we only show the main differences and a small part of the proof. First of all, we decompose the error: uk r−u(tk) = uk r−PRu(tk) + PRu(tk)−u(tk).(4.31) Now, we define ϑk:= uk r−PRu(tk) and %k:= PRu(tk)−u(tk),so: uk r−u(tk) = ϑk+%k.(4.32) If we use the triangle inequality, we have 1 Mkuk r−u(tk)k2≤2 M M X k=1 kϑkk2+2 M M X k=1 k%kk2.(4.33) The second term on the right hand side of (4.33) can be bounded by the properties of the projection (4.18)-(4.19) and the Poincar`e’s inequality (4.12): 1 M M X k=1 k%kk2≤C8 d X j=r+1 λj.(4.34) For the first term, we need to introduce some notation: for k= 1, . . . , m ¯ ∂ϑk:= ϑk−ϑk−1 ∆t, vk:= u(tk)−¯ ∂PRu(tk) = wk+zk, wk:= ut(tk)−¯ ∂u(tk), zk:= ¯ ∂u(tk)−¯ ∂PRu(tk). For all ψ∈Xr,we have (¯ ∂ϑk, ψ) + a(ϑk, ψ)=(¯ ∂uk r, ψ) + a(uk r, ψ)−(¯ ∂PRu(tk), ψ)−a(PRu(tk), ψ) = (f(tk), ψ)−(F(uk r), ψ)−(G(uk r), ψ)−(¯ ∂PRu(tk), ψ)−a(u(tk), ψ) = (vk, ψ)+(F(u(tk)) −F(uk r), ψ)+(G(u(tk)) −G(uk r), ψ). (4.35) Letting ψ:= ϑk,we get kϑkk2−(ϑk, ϑk−1)+∆tνk(ϑk)xk2= ∆t(vk, ϑk) + ∆t(F(u(tk)) −F(uk r), ϑk)+∆t(G(u(tk)) −G(uk r), ϑk).(4.36) We now have to estimate the r.h.s. of (4.36). The first one, if we use the Cauchy- Schwarz inequality, we obtain (vk, ϑk)≤ kvkkkϑkk.(4.37) 40 CHAPTER 4. THE POD FOR THE BURGERS EQUATION The second one, if we do the same as in the page 139 of [19] we obtain the following estimate: |(F(u(tk)) −F(uk r), ϑk)| ≤ νkϑkk2 H1+C9kϑkk2 L2+C10k%k2 H1.(4.38) For the last term, since it is not included in [19], we are going to use the same ideas as in [32]. We start by: (G(u(tk)) −G(uk r), ϑk) = (G(u(tk)) −G(PRu(tk)), ϑk) + (G(PRu(tk)) −G(uk r), ϑk).(4.39) Now, if we move the last term on the right hand side of (4.39) to the left hand side of (4.36) and apply the strong monotonicity property of G(4.15): (G(uk r)−G(PRu(tk)), ϑk)≥C11k(ϑk)xk3 L3.(4.40) If we apply the Lipschitz continuity of G(4.16) to the first term on the right hand side of (4.36): (G(u(tk)) −G(PRu(tk)), ϑk)≤C˜ Mk(%k)xkL3k(ϑk)xkL3,(4.41) where ˜ M= max{k(u(tk))xkL3,kPR(u(tk))xkL3}and we have used the Lemma (8). Using the first estimate in Lemma (6) and Young’s inequality: k(%k)xkL3k(ϑk)xkL3≤C2 3k(%k)xkk(ϑk)xkL2≤ ≤νk(ϑk)xk2 L2 2+1 2νC4 3k(%k)xk2 L2.(4.42) Using (4.36)-(4.42), we obtain the following approximation kϑkk2 L2−(ϑk, ϑk−1)+∆tνk(ϑk)xk3 L3≤ ≤ kvkkL2kϑkkL2+ (ν+C13 +ν 2)k(ϑk)xk2 L2+ (C10 +C4 3 2ν)k(%k)xk2 L2.(4.43) If we follow the rest of the proof in the page 140 of [19], choosing a small ∆tand summing over k, we obtain k(ϑk)xk2 L2≤eC14T(k(v0)xk2 L2+ ∆t k X j=1 (kϑk2 L2+C15k(%)xk2 L2)).(4.44) Finally, using the estimate of vkas done in [19] and (4.34) we obtain the main result (4.30). 4.3. PRACTICAL IMPLEMENTATION AND NUMERICAL RESULTS 41 4.3 Practical Implementation and Numerical Results The main goal of this section is to show the improvements of the AV-POD-G-ROM over the POD-G-ROM in terms of accuracy. We take the same initial condition and similar parameters used in [32, 21, 17]. In order to perform the computations, we consider the initial condition of (4.1) as : u0(x) =      1 if x∈(0,1 2], 0 if x∈(1 2,1). (4.45) The following numerical discretization has been used in all our numerical tests. The computational domain is Ω = (0,1) and the final time is T= 1.The parameters that we are going to use are: mesh size h=1 1024 for the Finite Element (FE) full order discretization; time step ∆t= 10−3; diffusion parameter ν= 10−5; number of snapshots M= 1001. Note that the used diffusion parameter is rather sharp with respect to the one used in [32, 21, 17], for which ν≥10−3. Firstly, we run a Direct Numerical Simulation (DNS). In addition, given the spatial and temporal resolution used, we can consider the resolution of the problem by finite elements as our benchmark solution, see Figure (4.2). The 1001 snapshots are generated by piecewise linear finite elements in space and backward Euler method in time, and they are used to generate the POD modes (we can see some of them in Figure (4.1)). The first 30 modes are used to generate the POD-G-ROM (4.4) and the AV-POD-G-ROM (4.7) with c= 10−4. In order to employ the DEIM method to the non-linear term, we consider 1001 snapshots too of it. We have to take a determinated number of mode in order to get more than the 99% of energy for it. In this case, if we consider ˜r= 60 we get it, but to show the efficiency of the DEIM algorithm we considered smaller values of ˜rtoo. We show the space-time evolution of the DNS in Figure (4.2), POD-G-ROM model in Figure (4.3) and AV-POD-G-ROM models with different values of ˜rin Figures 4.4-4.5-4.6. 48 CHAPTER 4. THE POD FOR THE BURGERS EQUATION Chapter 5 The POD for the Navier-Stokes equations In this chapter, we are going to study the numerical approximation of incompressible flows with the POD-ROMs. For this purpose, we consider the Navier-Stokes equations: ut−ν∆u+ (u·∇)u+∇p=fin Ω ×(0, T), ∇·u= 0 in Ω ×(0, T),(5.1) in a boundary domain Ω ⊂Rd, d ∈ {2,3},and with initial condition u(0) = u0. In the above equation, uis the velocity field, pthe pressure, ν > 0 the kinematic viscosity coefficient and fthe forcing term. It is necessary to complete (5.1) with boundary conditions, so for simplicity we are going to consider Dirichlet homogeneous condition. In order to give a variational formulation of problem (5.1), let us consider the velocity space: X= (H1 0)d={v∈(H1(Ω))d:v= 0 on ∂Ω}, and the pressure space: Q=L2 0(Ω) = {q∈L2(Ω) : ZΩ qdx = 0}. The weak formulation of (5.1) reads: 49 50 CHAPTER 5. THE POD FOR THE NAVIER-STOKES EQUATIONS Given f∈L2(H−1),find u: (0, T)→X, p : (0, T)→Qsuch that:      d dt(u, v)+(u·∇u, v) + ν(∇u, ∇v)−(p, ∇·v) = hf, vi ∀v∈Xin D0(0, T) (∇·u, q) = 0 ∀q∈Q u(0) = u0, (5.2) where h·,·i it is the duality between Xand its dual X0and D0(0, T) is the space of distributions in (0, T). Once obtained the weak formulation, we are going to give a FE approximation of (5.2). Let {Th}h>0be a family of regular triangulations of ¯ Ω.For any mesh cell k∈Th,the diameter will be called hkand h= maxk∈Thhk.The FE space for velocity and pressure are Xh⊂Xand Qh⊂Qrespectively. Therefore the FE approximation of (5.2) is: Find (uh, ph) : (0, T)→Xh×Qhsuch that:      d dt(uh, vh)+(uh·∇uh, vh) + ν(∇uh,∇vh)−(ph,∇·vh) = hf, vhiin D0(0, T) (∇·uh, qh) = 0 a.e. in (0, T) uh(0) = u0h, (5.3) for any (vh, qh)∈Xh×Qh,and uoh is some stable approximation to u0into Xh. Now, if we collect snapshots of the velocity field at certain time instances and we apply a standard POD-Galerkin ROM, we obtain an approximation of the solution of problem (5.1): u(x, t)≈ur(x, t) = r X j=1 aj(t)ϕj(x), where Xr=span{ϕ1, . . . , ϕr}is the space for the POD setting and {ai(t)}r i=1 are the time coefficients. Then, if we replace uwith urin (5.1), using the Galerkin method and projecting the resulted equations onto Xr,we obtain the POD-G- ROM for the Navier-Stokes equation (5.1): (∂ur ∂t , ϕ)+(ur·∇ur, ϕ) + ν(∇ur,∇ϕ) = 0 ∀ϕ∈Xr.(5.4) We note that, the POD modes are linear combinations of the snapshots, and thus belong to the space of discrete divergence-free functions. Also, without loss of generality, f≡0. 51 Remark 3. Since the POD velocity modes are linear combinations of the velocity snapshots, the POD velocity modes satisfy the boundary conditions in (5.1). This is because of the particular choice we have made at the beginning to work with homogeneous Dirichlet boundary conditions. In general, one has to manipulate the velocity snapshots set. This is the case, for instance, of steady-state nonhomogeneous Dirichlet boundary conditions, for which is preferable to consider a proper lift in order to generate POD velocity modes for the lifted velocity snapshots, satisfying homogeneous Dirichlet boundary conditions. This would lead to work with centered-trajectory method in the POD-ROM setting [16]. Remark 4. In (5.4), the pressure term vanishes due to the fact that the POD modes are weakly divergence-free. However, there exists different methods to recover the pressure, two of the most common methods are explained in [28]. These methods are used to solve the stability issues too, due to the fulfillment of the infsup condition in the reduced order framework. The first one enriches the velocity space adding more modes, which is called supremizer enrichment. The second one consists of adding a Possion equation for pressure and new boundary conditions for pressure too. We also consider another method, the Local Projection Stabilization(LPS) [24, 22]. This method leads to a coupled velocity-pressure POD-ROM that uses pressure modes as well to compute the reduced order pressure by adding a LPS term to the reduced equations (see Section 6.2.2). From the POD-G-ROM (5.4), we obtain an autonomous system for the coefficient a(t) : · a(t) = Aa +aCaT,(5.5) where A,C correspond to the linear and quadratic term in the numerical discretization of the Navier-Stokes equations, respectively. The finite dimensional system (5.5) can be written componentwise as follow: · ai(t) = r X j=1 Aijaj(t) + r X j=1 r X m=1 Cijmam(t)aj(t),(5.6) where Aij =ν(∇ϕi,∇ϕj),(5.7) Cijm =−(ϕj·∇ϕm, ϕi).(5.8) 52 CHAPTER 5. THE POD FOR THE NAVIER-STOKES EQUATIONS 5.1 S-POD-G-ROM The POD-G-ROM is an efficient tool for many applications, for example it works well on diffusion-dominated and laminar fluid flows. However, the POD-G-ROM for convection-dominated or turbulent flows yields inaccurate result, due to the effect of the discarded POD modes {ϕr+1, . . . , ϕd}on the retained ones, which is not taken into account. Indeed, although the disregarded modes do not contain a significant amount of the system’s kinetic energy, they have a significant role in the dynamics of the reduced-order system. For this reason, we are going to introduce a closure model to address this issue, the Smagorinsky POD-G-ROM (S-POD-G-ROM), see [25]. This closure model consists in adding, both to the offline and online stage, an eddy viscosity: νS:= (Csδ)2k∇ukF,(5.9) where Csis the Smagorinsky constant, δis the lengthscale and k∇ukFis the Frobenius norm of ∇u. Adding this into (5.3), we obtain the following problem:              Find (uh, ph) : (0, T)→Xh×Qhsuch that d dt(uh, vh)+(uh·∇uh, vh) + ν(∇uh,∇vh) +(Csδ)2(k∇uhkF∇uh,∇vh)−(ph,∇·vh) = hf, vhiin D0(0, T) (∇·uh, qh) = 0 a.e. in (0, T) uh(0) = u0h. (5.10) Now if we apply the same process for the online stage as in the previous section, we obtain the new autonomous system (S-POD-G-ROM): · a(t) = (A+˜ A)a+aCaT.(5.11) This finite dimensional system (5.11) can be written componentwise as follows: · ai(t) = r X j=1 (Aij +˜ Aij)aj(t) + r X j=1 r X m=1 Cijmam(t)aj(t)∀i= 1, . . . , r, (5.12) where ˜ Aij =−(Csδ)2(k∇urkF∇ϕj, ϕi).(5.13) As presented in [32], the S-POD-G-ROM has a disadvantage, which is to compute the matrix ˜ A(depending on k∇urkF) at every time step. Therefore, we 5.2. NUMERICAL STUDIES 53 propose to employ the DEIM algorithm (3.5.1) to reduce the computational complexity. So, we obtain: k∇ukF= ˜r X j=1 αj(t) ˜ϕj,(5.14) where ˜ X˜r= span{˜ϕ1, . . . ˜ϕ˜r}and αj(t) are the time coefficients. Consequently, if we replace (5.14) into our autonomous system (5.11), we obtain: · a=Aa +αTDa +aTCa. (5.15) This finite dimensional system (5.15) can be written componentwise as follow: · ai(t) = r X j=1 Aijaj(t)+ ˜r X j=1 r X m=1 Dijmαj(t)am(t)+ r X j=1 r X m=1 Cijmaj(t)am(t)∀i= 1, . . . , r, (5.16) where Dijm =−(Csδ)2( ˜ϕj∇ϕm,∇ϕi).(5.17) 5.2 Numerical Studies In this section, we present numerical results for the S-POD-G-ROM we introduced in the previous section. The numerical experiments are performed on the benchmark problem of the 2D unsteady flow around a cylinder with circular cross-section [14], and we take the same parameters as in [2] and [22]. 5.2.1 Setup for numerical simulations Following [29], the computational domain is given by a rectangular channel with a circular hole (see Figure (5.1)): Ω = {(0,2.2) ×(0,0.41)}\{x: (x−(0.2,0.2))2≤0.052}. 54 CHAPTER 5. THE POD FOR THE NAVIER-STOKES EQUATIONS Figure 5.1: Computational grid No slip boundary conditions are prescribed on the horizontal walls and on the cylinder, and a parabolic inflow profile is provided at the inlet u(0, y, t) = (4Umy(A−y)/A2,0)T, with Um=u(0, A/2, t)=1.5m/s and A= 0.41 mthe channel height. At the outlet, we impose outflow (do nothing) boundary conditions (ν∇u−pId)n= 0 with nthe outward normal to the domain. The chosen Reynolds number is Re = 100, depending on the inflow velocity U= 2Um/3 = 1m/s, the cylinder diameter D= 0.1m, and the kinematic viscosity ν= 10−3m2/s and there is no external forcing, f= 0 m/s2.In the fully developed periodic regime, a vortex shedding can be observed behind the obstacle, resulting in the well-known von K´arm´an vortex street, see Figure 5.2. For the evaluation of computational results, we are interested in studying the temporal evolution of the following quantities of interest. The first one is the kinetic energy of the flow, given by: EKin =1 2kuk2 L2. The others quantities of interest are the drag and lift coefficients. In the present work, they are computed as volume integrals [18]: cD=−2 DU [(∂tu, vD)+(u·∇u, vD) + ν(∇u, ∇vD)−(p, ∇·vD)], cL=−2 DU [(∂tu, vL)+(u·∇u, vL) + ν(∇u, ∇vL)−(p, ∇·vL)], for arbitrary test functions vD, vL∈H1such that vD= (1,0)Ton the boundary of the cylinder and vanishes on the other boundaries, vL= (0,1)Ton the boundary of the cylinder and vanishes on the other boundaries. In the actual computations we have used the approach in ([33]), where the pressure term is not necessary to compute cD, cL,since we used discretely divergence-free test functions vD, vL obtained by Stokes projection. 5.2. NUMERICAL STUDIES 55 Figure 5.2: Final velocity FOM. 5.2.2 FOM and POD modes Our aim is to compare the performance of the S-POD-G-ROM with the POD-G- ROM. The numerical method to get the snapshots is the Smagorinsky method (5.1), with a spatial discretization using P2−P1FE for the pair velocity-pressure on a relatively coarse computational grid (see Figure (5.1), h= 6,6×10−2), resulting in 6175 d.o.f for velocities and in 1735 d.o.f for pressure. For the time discretization, a semi-implicit Backward Differentation Formula of order 2 (BDF2) has been applied, which guarantees a good balance between numerical accuracy and computational complexity [1]. In particular, the discrete time derivative has been approximated by the operator D2 t: D2 tun+1 h=3un+1 h−4un h+un−1 h 2∆tn≥1,(5.18) and we have considered an extrapolation for the convection velocity: ˜un h= 2un h−un−1 h, n ≥1,in order to guarantee a second order accuracy in time. At the inizialization (n= 0), we have to consider u−1 h=u0 h=u0h,being u0hthe discrete initial condition, so the time scheme for the first iteration reduces from BDF2 to the semi-implicit Euler method for the first time step (∆t)0= (2/3)∆t. 56 CHAPTER 5. THE POD FOR THE NAVIER-STOKES EQUATIONS In the FOM simulation, the initial condition is a zero velocity field and the time step is ∆t= 2 ·10−3s. Time integration is performed till a final time T= 7s. For Re = 100, in the time period [0,5]sthe flow is expected to develop to full extent, including a subsequent relaxation time. We have plotted the time evolution of drag and lift coefficients in Figure (5.3) and kinetic energy in Figure (5.4) to show this. Figure 5.3: Temporal evolution of drag and lift coefficients. The POD modes are generated in L2by the method of snapshots with centeredtrajectories by storing every FOM solution from T= 5sto one period. The full period length of the statistically steady state is 0.332sfor Re = 100, so we collect 167 snapshots. In Figure (5.5) we can see the decay of the eigenvalues associated to the correlation matrix of the velocity, and for the DEIM algorithm applied to 5.2. NUMERICAL STUDIES 57 Figure 5.4: Temporal evolution of kinetic energy. Figure 5.5: Decay of the eigenvalues and captured energy. the eddy viscosity, and the captured energy calculated by 100 Pr k=1 λk/Pd k=1 λk where λkare the corresponding eigenvalues. Note that with 5 POD velocity modes, we already capture more than 99% of the energy, while for the DEIM algorithm we need around 20 modes. Once the POD modes are generated, the S-POD-G-ROM and POD-G-ROM are constructed as it is explained in section 5.1, using the same time discretization as the FOM and run in [5,7]sfor Re = 100 using only up to 7 POD velocity modes and 30 DEIM modes in order to get good results. We also decided to employ an adaptive in time algorithm for the Smagorinsky coefficient Cs, by adjusting its dissipation with the FOM energy, see next section. This caused a considerable improvement in the results. 64 CHAPTER 6. CONCLUSION AND OUTLOOK 6.2.2 Pressure recovery and turbulent flows In chapter 5, we already mentioned several methods to recover the pressure in fluid dynamics reduced order simulations. A recent proposal has been to apply a Local Projection Stabilization (LPS) technique in the POD-ROM framework [24, 22]. This method is a velocity-pressure ROM that uses pressure modes as well to compute the reduced order pressure needed for instance in the computation of relevant engineering quantities. It circumvents the standard discrete inf-sup condition for the POD velocity-pressure spaces, whose fulfillment can be rather expensive in realistic applications in Computational Fluid Dynamics (CFD). Also, the velocity modes do not have to be either strongly or weakly divergence-free, which allows the use of snapshots generated, for instance, with penalty or projection-based stabilized methods. In particular, in my ongoing research (see Appendix A), I obtained some new results combining the LPS for pressure stability and the RBF technique (see section 3.5.2) to approximate the eddy viscosity, to tackle the increase of Reynolds number in multi-parametric POD-ROMs [13]. We believe this is a first step towards the challenging simulation of realistic turbulent flows by ROMs, which would be the goal of future research. Appendix A Ongoing research Figure A.1: Mesh of the problem (Flow past cylinder). Figure A.2: Decay of POD eigenvalues of velocity, pressure and eddy viscosity. 65 66 APPENDIX A. ONGOING RESEARCH Figure A.3: Kinetic energy relative error for different Reynolds values including Re = 187.5 (not considered in the snapshots sample). Figure A.4: Temporal evolution of L2norm of relative error of velocity (right) and pressure (left). Figure A.5: Drag and Lift coefficient for Re = 187.5. 67 Figure A.6: Comparison of the velocity field for FOM (left) and ROM (right) for Re = 187.5.The fields are depicted for different time instant equal to t = 7.5 s , 8 s , 8.5 s and 9s and increasing from top to bottom. The ROM solutions are obtained with 8 modes for velocity and pressure and 60 modes for eddy viscosity. 68 APPENDIX A. ONGOING RESEARCH Figure A.7: Comparison of the pressure field for FOM (left) and ROM (rigth) for Re = 187.5.The fields are depicted for different time instant equal to t = 7.5 s , 8 s , 8.5 s and 9s and increasing from top to bottom. The ROM solutions are obtained with 8 modes for velocity and pressure and 60 modes for eddy viscosity. Bibliography [1] N. Ahmed and S. Rubino,Numerical comparisons of finite element stabilized methods for a 2d vortex dynamics simulation at high reynolds number, Computer Methods in Applied Mechanics and Engineering, (2019), pp. 191– 212. [2] B. Archilla, J. Novo, and S. Rubino,Error analysis of proper orthogonal decomposition data assimilation schemes for the navier-stokes equations, (2020). [3] M. Aza¨ ıez, T. C. Rebollo, and S. Rubino,Streamline derivative projection-based pod-rom for convection-dominated flows. part i : Numerical analysis, 2017. [4] M. Barrault, Y. Maday, N. Nguyen, and A. Patera,An ‘empirical interpolation’ method: Application to efficient reduced-basis discretization of partial differential equations, Comptes Rendus Mathematique, 339 (2004), p. 667–672. [5] M. C. Bramwell,Introductory functional analysis with applications, The Mathematical Gazette, 63 (1979), p. 137–138. [6] S. Chaturantabut and D. C. Sorensen,Discrete empirical interpolation for nonlinear model reduction, in Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference, 2009, pp. 4316–4321. [7] L. Cordier,R´eduction de dynamique par d´ecomposition orthogonale aux valeurs propres ( pod 1 ), 2006. [8] J. Edmonds,Matroids and the greedy algorithm, Mathematical programming, 1 (1971), pp. 127–136. [9] L. Evans,Partial Differential Equations, Graduate studies in mathematics, American Mathematical Society, 2010. 69 70 BIBLIOGRAPHY [10] K. P. F.R.S.,Liii. on lines and planes of closest fit to systems of points in space, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 2 (1901), pp. 559–572. [11] G. Golub and C. Van Loan,Matrix Computations, Johns Hopkins Studies in the Mathematical Sciences, Johns Hopkins University Press, 2013. [12] N. J. Higham,Matrix nearness problems and applications, 1989. [13] S. Hijazi, G. Stabile, A. Mola, and G. Rozza,Data-driven pod-galerkin reduced order model for turbulent flows, Journal of Computational Physics, 416 (2020), p. 109513. [14] E. H. Hirschel,Flow simulation with high-performance computers ii, AIAA Journal, 35 (1996), pp. 1560–1560. [15] C. Homescu, L. R. Petzold, and R. Serban,Error estimation for reduced-order models of dynamical systems, SIAM Rev., 49 (2007), p. 277–299. [16] T. Iliescu, V. John, S. Schyschlowa, and D. Wells,Supg reduced order models for convection-dominated convection-diffusion-reaction equations, 09 2014. [17] T. Iliescu and Z. Wang,Are the snapshot difference quotients needed in the proper orthogonal decomposition?, 2013. [18] V. John,Reference value for drag and lift of a two-dimensional time dependent flow around cylinder, International Journal for Numerical Methods in Fluids, 44 (2004), pp. 777 – 788. [19] K. Kunisch and S. Volkwein,Galerkin proper orthogonal decomposition methods for parabolic problems, Numerische Mathematik, 90 (2001), pp. 117– 148. [20] Z. Luo, J. Chen, I. M. Navon, and X. Yang,Mixed finite element formulation and error estimates based on proper orthogonal decomposition for the nonstationary navier–stokes equations. [21] C. Mou, B. Koc, O. San, L. G. Rebholz, and T. Iliescu,Data-driven variational multiscale reduced order models, Computer Methods in Applied Mechanics and Engineering, 373 (2021), p. 113470. [22] J. Novo and S. Rubino,Error analysis of proper orthogonal decomposition stabilized methods for incompressible flows, SIAM Journal on Numerical Analysis, 59 (2020), pp. 334–369. BIBLIOGRAPHY 71 [23] A. Quarteroni, A. Manzoni, and F. Negri,Reduced basis methods for partial differential equations: An introduction, 01 2015. [24] S. Rubino,Numerical analysis of a projection-based stabilized pod-rom for incompressible flows, 2019. [25] P. Sagaut and Y.-T. Lee,Large eddy simulation for incompressible flows: An introduction. scientific computation series, Applied Mechanics Reviews, 55 (2002), pp. 115–. [26] L. Sirovich,Turbulence and the dynamics of coherent structures. i - coherent structures. ii - symmetries and transformations. iii - dynamics and scaling, Quarterly of Applied Mathematics - QUART APPL MATH, 45 (1987). [27] L. SIROVICH,Turbulence and the dynamics of coherent structures part iii: Dynamics and scaling, Quarterly of Applied Mathematics, 45 (1987), pp. 583– 590. [28] G. Stabile and G. Rozza,Finite volume pod-galerkin stabilised reduced order methods for the parametrised incompressible navier–stokes equations, Computers and Fluids, 173 (2018), pp. 273–284. [29] S. Turek,Recent benchmark computations of laminar flow around a cylinder, (2000). [30] S. Volkwein,Proper orthogonal decomposition: Theory and reduced-order modelling, Lecture Notes, University of Konstanz, (2012). [31] , Proper orthogonal decomposition methods for partial differential equations, SIAM Review, 63 (2021), pp. 242–244. [32] Z. Wang,Reduced-Order Modeling of Complex Engineering and Geophysical Flows: Analysis and Computations, PhD thesis, Faculty of the Virginia Polytechnic Institute and State University, 4 2012. [33] C. Zerfas, L. Rebholz, M. Schneier, and T. Iliescu,Continuous data assimilation reduced order models of fluid flow, Computer Methods in Applied Mechanics and Engineering, 357 (2019), p. 112596.