Reduced Basis Method applied to the Smagorinsky Turbulence Model

Programa de Doctorado "Matemáticas" PHD DISSERTATION Reduced Basis Method applied to the Smagorinsky Turbulence Model Author Cristina Caravaca García Supervisors Prof. Dr. Tomás Chacón Rebollo Prof. Dra. Macarena Gómez Mármol Sevilla, 2022 ii En memoria de mi abuelo Pedro. A mis padres y a mi hermano, siempre luchando. A todos los que creyeron en mi. iv Acknowledgements Primero, agradecer a mis directores Macarena y Tomás por darme la oportunidad de vivir esta experiencia. Sobre todo por estar conmigo hasta el final, por la infinita paciencia que han tenido conmigo y el apoyo en los momentos no tan buenos. A mis compañeros de despacho, que han sido muchos. En especial a Enrique por saber ayudarme en mis inicios y a Javi por estar ahí siempre que lo he necesitado. También mencionar a los compañeros del IMUS que me han acogido como una más. Al departamento de EDAN y al grupo de Investigación M2S2M. Muchos me dieron clase en mi etapa de estudiante de Grado y he tenido la gran suerte de trabajar y aprender con ellos. Incluso conocerles un poco mejor en congresos y en esas comidas de Navidad en Marzo. I acknowledge Prof. F. Hecht for receiving me at the LJLL, helping me to improve my research every day, and even teaching me some French. To Anouk, Elise, Jules and Gabriela, for helping me with my every day in the laboratory and for doing my stage a lot more pleasant. A mi familia, siempre comprensivos, ayudándome en lo que podían. En especial, a Ana y Antonio por acogerme en su casa como si fuera mía durante los últimos meses. Y por supuesto a mis padres, que siempre tuvieron para mí un “tu puedes”. Y junto a ellos, a mi hermano, por servirme de inspiración para luchar. Y por último, quiero agradecer a Javi su infinita paciencia y su apoyo incondicional. Sin él, no hubiese llegado hasta aquí. Resumen Esta tesis se enmarca en la resolución numérica de modelos que simulan el comportamiento de flujos turbulentos mediante técnicas de orden reducido. La resolución numérica de este tipo de flujos es compleja y costosa, a la par que necesaria para el diseño óptimo de edificios eco-eficientes, entre otras muchas aplicaciones en ingeniería y arquitectura. La incorporación de técnicas de Orden Reducido es clave para la reducción en órdenes de magnitud el tiempo y coste computacional en la resolución numérica de estos problemas. En esta tesis utilizaremos el modelo de Smagorinsky, un modelo de turbulencia de tipo LES (Large Eddy Simulation) basado en las ecuaciones de Navier-Stokes que permite la resolución de flujo turbulento sin la necesidad de mallas muy finas. Aún así, el coste computacional es elevado, sobre todo en casos 3D. Dentro de la modelización de orden reducido, existen varias técnicas que permiten la obtención del modelo reducido. En esta tesis nos centramos en el desarrollo de métodos de Bases Reducidas (BR). Usaremos de forma auxiliar la “Proper Orthogonal Decomposition” (POD). Para la obtención del modelo reducido, es necesario el cálculo del error entre el modelo aproximado y el modelo reducido mediante el método BR, lo cuál puede llegar a ser muy costoso y complejo. Es por ello que se estudia la deducción de estimadores a posteriori que permiten estimar este error y cuyo cálculo sea más rápido. El modelo de Smagorinsky es no lineal ya que deriva de las ecuaciones de Navier-Stokes y el desarrollo de estimadores para problemas no lineales requiere la utilización de técnicas matemáticas adaptadas. Para un modelo estacionario no lineal, encontramos estimadores basados en la teoría Brezzi-Rappaz-Raviart (BRR) de aproximación de ramas no singulares de problemas no lineales paramétricos. Es una teoría esencialmente basada sobre el teorema de la función implícita. Este estimador ha sido ya elaborado para flujos estacionarios, aunque solo se ha aplicado a flujos 2D. Por lo que, comenzamos aplicando este estimador a un caso 3D, obteniendo grandes reducciones en cuanto a tiempo y coste computacional. También aplicamos este estimador a un problema realista de diseño de claustros orientado a la optimización del confort térmico en la planta baja. Los parámetros a considerar para el problema son la altura y anchura del pasillo que rodea el claustro. Conseguimos obtener una respuesta óptima analizando 625 posibilidades en 16 minutos. Uno de los principales desafíos que abordamos es extender la obtención de estimadores a posteriori a problemas no estacionarios. Para empezar, es necesario realizar estudios previos sobre estimaciones a priori que involucren a la velocidad y presión. Así, somos capaces de desarrollar un estimador a posteriori asumiendo ciertas hipótesis. Por otro lado, desarrollamos una alternativa basándonos en la teoría de turbulencia en equilibrio estadístico de Kolmogórov, por la cuál sabemos que existe una cascada de energía que se propaga desde los grandes torbellinos hacia los más pequeños disipando la energía en ellos gracias a la viscosidad. Esta cascada genera un espectro de energía inercial con una forma determinada que utilizamos como estimador. Validamos dicho estimador con un test académico, utilizando una estrategia POD+Greedy, obteniendo resultados similares utilizando el estimador y el error real cometido. Abstract This PhD dissertation addresses the numerical simulation of models that simulate the behavior of turbulent flows through reduced-order techniques. The numerical simulation of this kind of flows is complex and pricey, as well as necessary to the Eco-efficient building optimized design, among many others applications in engineering and architecture. The integration of reduced-order techniques is the key to reducing by orders of magnitude the time and computational cost in the numerical simulation of these problems. In this dissertation, we use the Smagorinsky model, a Large Eddy Simulation (LES) model based upon the Navier-Stokes equations that allows for the resolution of turbulent flows with coarser meshes. Even then, the computational cost is high, especially in 3D cases. With respect to the Reduced Order Model (ROM), there exist some techniques to obtain the ROM. In this dissertation, we focus on the development of Reduced Basis (RB) method. Upon an ancillary basis, we shall use Proper Orthogonal Decomposition (POD). For ROM obtainment, it is necessary to compute the error between the approximated model and the reduced model through the RB method, which could become very expensive and complex. This is the reason to study a posteriori estimates to estimate the error and which computation is faster. The Smagorinsky model is non-linear since it is derived by the Navier-Stokes equations and the estimator development for non-linear problems requires the use of adapted mathematical techniques. For steady non-linear models, we find estimates based on the Brezzi-Rappaz-Raviart (BRR) theory of non-singular branches approximation of parametric non-linear problems. It is a theory essentially based upon the Implicit function theorem. This estimator has already been developed for steady flows, although it has only be applied to 2D flows. Therefore, we shall start applying this estimator to the 3D case, obtaining large reductions in time and computational cost. We also apply this estimator to a realistic design problem for a cloister focused on the thermal comfort optimization of the ground floor. The parameters to consider for the problem are the height and width of the corridors around the cloister. We obtain an optimal answer by analyzing 625 possibilities in 16 minutes. xvi Notation and Abbreviations PDE ...... Partial Differential Equation. NS ....... Navier-Stokes. LES ....... Large Eddy Simulation. FE ....... Finite Element. EIM ...... Empirical Interpolation Method. RBF ...... Radial Basis Function. ROM ...... Reduced-Order Model. RB ....... Reduced Basis. POD ...... Proper Orthogonal Decomposition. BRR ...... Brezzi-Rappaz-Raviart. Variable Name SI units ρmass density kg/m3 µdynamical viscosity kg/ms νkinematic viscosity m2/s κthermal diffusivity m2/s cpspecific heat capacity J/kgK kthermal conductivity W/mK αconvective heat transfer W/m2K βthermal expansion (1/θfor ideal gases) K−1 u,v,z,wvelocity field m/s E=1 2|u|2kinetic energy per mass density m2/s2 ε=2ν|Du|2dissipation per mass density m2/s3 ∇udeformation tensor s−1 Du= (1/2)(∇u+∇ut)symmetric deformation tensor s−1 fsource term per mass unity m/s2 p,q,rpressure per mass density m2/s2 Ucharacteristic velocity m/s Lcharacteristic length m Table 1: Table of coefficients and variables. Introduction Nowadays, many of the problems in the fields of physics and engineering can be characterized by Partial Differential Equations (PDEs). At the same time, PDEs involve parameters that change from a resolution to another, modifying the results. This derives in the resolution of highly complex parametric systems, such as evolutionary and inverse problems, data assimilation in non-linear PDEs and optimal design of complex devices, among others. Along this dissertation, we focus on Computational Fluid Dynamics, which is a branch of fluid mechanics. These flows frequently depend on parameters such that the ratio between high and width of a building to study its aerodynamics, the attack angle of an airfoil, or governing physical parameters such that viscosity, density, thermal conductivity, etc. PDEs are solved applying mathematical methods of numerical approximation. Design problems or fluid variability calls for the solution of PDEs on many occasions, changing the parameters, what can become very computationally expense. Even the use of the use of High Performance Computation (HPC) could take a considerable time. Hence, the need of techniques to develop Reduced Order Models (ROM) whose resolution is expected to be much faster and with less resources. The main purpose of the ROM method is to reduce by orders of magnitude the computational cost of numerical simulation of parameterized partial differential equations. This method is able to catch the dominant characteristics of the problem accepting some admissible errors. This method drastically reduces the calculation and allows to carry out realistic simulations in acceptable computing time without the need of HPC. The ROMs are constructed in an off-line phase, which requires solving the full-order model (FOM) for several well-chosen cases. This information is treated to capture the most important one, using appropriate techniques, and to build the ROM that is run with highly reduced computational time in an on-line phase. Some basic references about ROM are the works of Quarteroni et al. [ 38 ], Schilders et al. [ 45 ] and Chinesta et al. [ 13 ]. There exists various kinds of ROMs to solve parametric PDEs. In this dissertation, we shall mainly develop Reduced Basis (RB) methods to compute parametric turbulent 1 2Introduction flows. We shall also use the Proper Orthogonal Decomposition (POD) as an auxiliary ROMs. These techniques will be introduced in Chapter 1. The use of RB method requires the design of a posteriori error bound estimators, which is the main purpose of this dissertation. The RB method is built by means of the greedy algorithm, which looks for an optimal basis function to add to the current RB space. The estimator usually is based upon the dual norm of the residual that allows certifying the method. For linear PDEs, the procedure to obtain an a posteriori error bound estimator is well-known, applying continuity and coercivity results to the error (between full order and reduced order solutions) equation (cf. [ 38 ]. For non-linear PDEs, the a posteriori error bound estimator can be developed following the Brezzi-Rappaz-Raviart (BRR) theory (cf. [ 8 ]) for approximation of regular branches of solutions of non-linear PDEs. The POD is the most used technique for the tensorized approximation of bi-parametric functions, in this case, the full order solution of evolution PDEs. This technique looks for minimizing the mean square error between the full order solution and the reduced solution over all possible subspace of the “snapshots” data of given dimension. This space is formed by the orthogonal functions (POD basis) that are obtained by singular value decomposition (SVD) of a suitable correlation matrix. Once the reduced space is built, a Galerkin projection of the governing PDEs can be employed to obtain a low-order dynamical system, in to the parameter window used to create the POD basis. The resulting low-order model is named standard POD-Galerkin ROM. As we mentioned before, we are interested in the application of ROM to turbulent flows. Most of the phenomena which carry on a difficulty in the numerical resolution are turbulence flows. This phenomenon is characterized by chaotic, multi-scale dynamics, both in space and time. The governing PDEs are the Navier-Stokes (NS) equations whose resolution becomes unaffordable for many real flows of interest. Using the mesh size and time steps to accurately solve a single flow past a 3D obstacle with well-developed turbulence would require dozens of HPC nodes (cf. [ 12 ]). To alleviate this, alternative models have been considered, such as Large Eddy Simulation (LES) models. The LES models are based upon the Kolmogórov’s theory of statistical equilibrium, introduced in 1941 by Andrey Kolmogórov. This theory describes how energy is transferred from large to smaller eddies, ( cf. [ 27 , 28 ]). This is at high Reynolds number (a physical parameter that triggers turbulence), the flow exhibits an energy cascade: large scale eddies are broken down into smaller and smaller eddies until the scales are fine enough so that viscous forces can dissipate their energy. LES model strategy is based upon to solve large grid scales and to approximate small scales Introduction 3 by an eddy diffusion term in the NS equations, which involves the mixing length of the sub-grid eddies, which is proportional to the grid size. Even the computation of this kind of models is quite complex. Most of the turbulence effects are genuinely 3D (mainly due to vortex stretching triggered by inertial effects) do not admit a simplification to 2D. The solution of parametric turbulent flows would call for the resolution of this model multiple times, which should be computed in affordable times. At this point, we dip into ROM for the LES model, in particular, for the Smagorinsky model. This model was introduced by Joseph Smagorinsky in 1963 (see [ 46 ]), adding an eddy viscosity term to the averaged Navier-Stokes equations. The primary objective of this dissertation is to build a posteriori error bounds for the finite element - reduced order solution of the Smagorinsky turbulence model. The use of BRR theory for the steady Smagorinsky model was developed by E. Delgado in his PhD dissertation [ 14 ]. The estimator has been validated for steady 2D problems only, therefore, as a first contact with this estimator, in this thesis, we validate the estimator with a numerical test in 3D considering the Reynolds number as the parameter for the problem. The results are shown in Section 1.5, in which we consider the lid-driven cavity problem for relatively low Reynolds numbers, for which the flow is steady. The outcomes are promising reducing the computational time from almost 2 hours to over 2 seconds, with speed-ups around 3000 . The remaining thesis is divided into two parts. In the first part, we seek for the application of a posteriori error bound estimator to a realistic problem of multiparameter design. In the second part, we seek to extend to the unsteady Smagorinsky model the techniques to construct RB solvers developed for steady problems. First, it is necessary to develop a priori estimates for the unsteady Smagorinsky model to prove well-possednes. Subsequently, we develop an a posteriori error estimator following the BRR theory. Lastly, we seek for an a posteriori estimate based on the statistical equilibrium turbulence theory developed by Andréi Kolmogórov in 1941, and we apply it to an academic test. Part I The use of courtyards in the Andalusian architecture goes back to the Roman and Muslim culture. They were used as a communal space, however, the courtyards have many others utilities. Nowadays, we can find this kind of structure in many buildings and their purpose goes from providing light into buildings to create a fresh environment in the center of the building to refrigerate the rest of it. The air movement will renew the 4Introduction air inside the courtyard, which will lower the temperature with respect to the outside and so the building temperature (cf. [ 16 , 32 ]). Cloisters are a kind of courtyard widespread in old buildings such as cathedrals and monasteries. In Figure 1, an example of this cloister is shown. We are interested in the optimal thermal comfort design of the cloister, studying the dynamics of the air that flows into the courtyard by forced convection model, as a first approximation to our problem. Figure 1: Merced Convent, placed in Córdoba (Spain). For the building design, it is interesting to know the distribution of temperature along the courtyard to select the best thermal comfort design. We consider two parameters, the width and the height of the corridor around the cloister. The resolution of this kind of flow could take a considerable computing time. This is the reason to consider to apply the Reduced Basis (RB) Method that we develop along Chapters 2 and 3. As a first approximation, we consider the fluid dynamics only in 2D. We develop the RB problem in Chapter 2 building an a posteriori error bound estimator following the BRR theory for velocity and pressure. In Chapter 3, we add the temperature to the problem using a forced convection model. We build the RB problem developing an a posteriori error bound estimator for temperature, following now standard techniques since the problem is linear. We are able to obtain the distribution of velocity, pressure, and temperature in barely 1.5 second in stead of 3.5 minutes with errors more than admissible. Introduction 5 Taking into account that our parameters are geometrical and the computational time for the ROM, it makes sense to apply our RB model to a design problem. We study which is the best combination of both geometrical parameters so that the temperature in the ground floor is the closest to the comfort one, in order to optimize the comfort properties of our cloister. We design a functional to minimize the temperature in the ground floor and we computed the value of this functional for 625 pairs of parameters in less than 16 minutes. We obtain an optimum design that optimizes the thermal comfort in the lower part (to be occupied by people) of the cloister. Part II The main goal of this part is to build RB based upon error estimators for the unsteady Smagorinsky turbulence model. We start this part studying a priori estimates for velocity and pressure for the time-space discrete Smagorinsky model in Chapter 4. We base this chapter on the results obtained by Volker John in Chapter 6 [ 25 ], where he developed a priori estimates for finite element solution of LES models. We obtain a priori estimates over appropriate space-time Sobolev spaces using inverse inequalities. This results are the key to apply the BRR theory in the next chapter. Secondly, in Chapter 5, we are able to apply the BRR theory for the development of a posteriori error bound estimator, not without finding difficulties. The key is the definition of the norm related to the problem. Intuition could say that we need a norm relating the velocity derivative in time, deformation and pressure. However, we can not ensure the inf-sup condition considering this norm. This is why we consider a norm related to velocity and pressure, we prove well-posedness depending on the time step chosen, the inf-sup condition is achieved, and we obtain the a posteriori error bound estimator. Unfortunately, the computation of the Stability Factor, necessary for the computation of the a posteriori estimator, is not low-cost, and it is necessary the implementation of further techniques to reduce the complexity. Finally, we look for an alternative, and we dip into the statistical equilibrium Kolmogórov theory. We explain this theory along Chapter 6, explaining the concept of energy spectrum and how to compute it. Lastly, we introduce the a posteriori error estimator using the energy spectrum. The continuous problem should achieve an specific energy spectrum, since the full-order model is intended to be a good approximation of the continuous 6Introduction problem, it should also achieve this energy spectrum in the resolved part of the inertial spectrum. The key is to use the error committed in this energy spectrum by the RB solution to select the new basis functions by the Greedy algorithm. To validate the estimate, we develop an academic test in which we compare the use of the Kolmogórov estimator with the use of the error between the full order and reduced order solution to develop the RB problem. The test is based on a problem for which we obtain a velocity field that satisfies the −5/3 law, and we build a RB problem using the POD+Greedy strategy. The number of bases and the error are similar to those obtained using the true error as “error estimator”, which indicates that the estimate works in this case. In this case, the speed-ups are around 20 , which is more than admissible as using the exact error does not provide a further speed-up improvement. FreeFem++ and post processing. All numerical test have been coded in FreeFem++ v. 4.8 (cf. [ 23 ]). Due to the complexity of the numerical test, we use parallel computation using: • PETSc package (Portable, Extensible Toolkit for Scientific Computation). This package is embedded in FreeFem++ and allows us to solve variational problems. For that purpose, it is necessary a domain partition depending on the processor number, which can be boarded by commands defined in the macro_ddm.idp file, included in FreeFem++. Then, we define the FE matrices and we solve the problem using the PETSc package, obtaining a local solution in each partition. Finally, we recover the solution defined in the complete domain. This process has been automated using macros. • MPI (Message Passing Interface). Again, we have commands that allow us the transmission of information between processors thanks to macro_ddm.idp file. We use these commands to compute independent loops in parallel, necessary for example, for the construction of the RB matrices. The main idea is to assign the computation of some loop iterations to each processor and finally share the result within the processors. We can apply this technique to the construction of matrices since these involve two loops with independent iterations. Introduction 7 Every offline phase along this dissertation has been performed in the cluster Anonimus2021 which allows a High Performance Computation (HPC). This cluster is composed by 18 nodes, two of them dedicated to computation on GPUs. The rest of the nodes are composed by CPUs AMD EPYC 7542. For the figures presented along this dissertation, we have used MATLAB 2016b (cf. [ 48 ]) and Paraview v. 5.9 (cf. [ 3 ]) 8 1 Basics of Reduced Basis modelling In this chapter, we introduce basic results that will be used along this dissertation. We introduce a general notation for Partial Differential Equations (PDEs) in Section 1.1, the Reduced Order Model (ROM) for PDEs in Section 1.2. We introduce the Empirical Interpolation Method (EIM) for the reduced linearization of some terms with respect to the parameter in Section 1.3. Finally, we show a numerical test, application of the Reduced Basis Method to the 3D lid-driven cavity flow. This test is an extension to the numerical test performed by E. Delgado et al. in [ 11 ]. In this paper, E. Delgado et al. develop a technique to construct a posteriori error estimator necessary for the Reduced Basis Method, and they validate the estimator for the 2D case. We apply this theory to the 3D case in Section 1.4. 1.1 Parametric Partial Differential Equations We shall consider in this dissertation parametric PDEs. Let us denote by µ µ µ= (µ1,...,µP) the parameter vector where P is the total number of parameters. We note that each µi for i=1,...,P can describe a property in the model, such as viscosity, density, thermal conductivity, etc. (physical parameters) or it can parameterize the shape of the domain (geometric parameters). Let us denote by D⊂RP the parametric set that will contain all possible parameter values. Finally, let us denote by Ω=Ω(µ µ µ)∈Rd , d=2,3 the spatial domain and If= (0,Tf) with Tf>0 the time domain, Γ the boundary of Ω and QT=If×Ω . 9 16 1.2. Model-Order Reduction Algorithm 4 POD+Greedy Set ε1,tol,ε2,tol >0, Nmax ∈N,µ µ µ∗∈Dtrain,Z= [ ] and S={ }; while N<Nmax do S=S∪{µ∗}; Compute uk h(µ∗)for k=1,...,L; Build S= [u1 h(µ µ µ∗),u2 h(µ µ µ∗),...,uL h(µ µ µ∗)]; [ξ1,...,ξM]=POD(S,Xh,ε1,tol); Z= [Z,ξ1,...,ξM]; [ϕ1,...,ϕN]=POD(Z,Xh,ε2,tol); XN={ϕi}N i=1; µ µ µ∗=argmaxµ µ µ∈Dtrain ∆N(µ µ µ,L); εN=∆N(µ µ µ∗,L); if εN≤εtol then Nmax =N; end if end while Algebraic Construction Let us recall that XN is the space generated by the RB functions {ϕi}N i=1 . The reduced solution uk N(µ µ µ)∈XN of (1.6) for any k=1,...,L can be described as a linear combination of {ϕi}N i=1 , uk N= N ∑ i=1 (uk N)i(µ µ µ)ϕi where uk N is the solution of the reduced linear system (MN+AN)uk N=fk N+MNuk−1 N,(1.7) being MN,AN∈MN×N(R) and fk N∈RN associated to the left and right hand side of (1.6) . These matrices depend on the parameter µ µ µ and their construction can take a long computing time depending on N since they are full matrices. To save that time in the online phase, we look for a linear representation of the problem (1.6) with respects to the parameters. This linear representation will allow us to define RB parameter-independent matrices that we store at the offline phase end. Remark 1.1. For the Greedy Algorithm 1-2, the construction of the spaces XN is recursive, this is, in order to build the N -iteration, we need the N−1 -iteration. The construction of 1. Basics of Reduced Basis modelling 17 the RB matrices could be done in a recursive way too. For example, for the matrix related to the M operator, given MN−1, then, MN=     (ϕN,ϕ1)Ω/∆t MN−1. . . (ϕN,ϕN−1)Ω/∆t (ϕ1,ϕN)Ω/∆t... (ϕN−1,ϕN)Ω/∆t(ϕN,ϕN)Ω/∆t      . Physical vs. geometrical parameters In terms of parameters, we shall identify two kinds, physical or geometric. Dimensionless numbers, as those described in Appendix A.4, are meant to be a physical parameter. In this case, the parameter representation is almost direct. However, it could be non-linear and linearization techniques should be applied. In the next section, we will introduce the EIM method that will help us to deal with this situation. For the geometric parameter, the dependency is intrinsic to the domain. In this case, we are able to fix a parameter µ µ µr∈D which defines a reference domain Ωr=Ω(µ µ µr) and we define problem (1.5) over Ωr instead of Ω(µ µ µ) thanks to a change of variable. We will see an application in Chapter 2, however, we refer to [ 41 , 42 , 38 , 14 ] for more details. 1.3 Empirical Interpolation Method The aim of the Empirical Interpolation Method (EIM) is to build a linear approximation gM of non-linear function g respect to the time t and paramter µ µ µ as follows gM(x;t,µ µ µ) = M ∑ i=1 σi(t,µ µ µ)qi(x),∀x∈Ω,t∈If,µ µ µ∈D for σi and qi for i=1,...,M to be determined. To build this approximation, let consider the space G={g(·;t,µ µ µ),µ µ µ∈D,t∈If} as the family of parameter-time-dependent functions. Then, the key of the EIM is to find a reduced approximation space GM=span{qi}M i=1⊂G in a first offline phase. The offline phase could be afforded using the Greedy algorithm. The coefficients {σi}M i=1 are obtained by the solution of M ∑ j=1 qj(xi)σj(t,µ µ µ) = g(xi;t,µ µ µ),i=1,...,M 18 1.4. ROM application to steady Smagorinsky turbulence model where TM={xi}M i=1 are the interpolation or magic points (cf. [ 33 ]). That is, in matrix form Bσ σ σ(t,µ µ µ) = gM(t,µ µ µ),∀t∈If,µ µ µ∈D(1.8) where (B)i j =qj(xi),(σ σ σ(t,µ µ µ))i=σi(t,µ µ µ),(gM)i=g(xi;t,µ µ µ), for i,j=1,...,M . The offline and online phases are described in Algorithm 5. Algorithm 5 EIM (computable version) function (B,T,G,S)=EIMOFFLINE(g,εEIM,Mmax,DEIM,∆t,L) Set S={},T={ },G={ }; (k1,µ µ µ1) = argmaxµ µ µ∈DEIM,k=0,...,L∥g(·;k∆t,µ µ µ)∥∞; r1=g(·;k1∆t;µ µ µ1); M=1; while M<Mmax do xM=argsupx∈Ω|rM(x)|; qM=rM/rM(xM); S=S∪{(kM,µ µ µM)}; T=T∪{xM}; G=G∪{qM}; Bi j =qj(xi),1≤i,j≤M; (kM+1,µ µ µM+1) = argmaxµ µ µ∈DEIM,k=0,...,L∥g(·;k∆t,µ µ µ)−gM(·;k∆t,µ µ µ)∥∞; rM+1=g(·;kM+1∆t,µ µ µM+1)−gM(·;kM+1∆t,µ µ µM+1),ε=∥rM+1∥∞; M=M+1; if ε≤εEIM then Mmax =M; end if end while end function function σ=EIMONLINE(B,t,µ µ µ,{x1,...,xM}) for i=1 : Mdo gi= (xi;t,µ µ µ); end for σ=B−1g; end function 1.4 ROM application to steady Smagorinsky turbulence model In [ 11 ], E. Delgado et al. have built a RB Smagorinsky model for the steady case. They introduce a posteriori error bound estimator for non-linear problems using the BRR 1. Basics of Reduced Basis modelling 19 theory (cf. [ 8 ]). In this section, we will introduce the basics of the theory developed by E. Delgado, which will be the foundation of Chapter 2 and 5. Let us consider the dimensionless Smagorinsky turbulence model. It is a Large Eddy Simulation (LES) model and it was first introduced by Joseph Smagorinsky in 1963 (see [ 46 ]). LES modelization is based upon the Kolmogorov Theory of statistical equilibrium theory that we will describe in detail in Chapter 6. In essence, we suppose that the Reynolds number is sufficiently high to ensure that energy is transferred from large to small scales, which are responsible for dissipation only. This generates an inertial spectrum in which only inertial forces are important, delimited by two scales which form the inertial range. The LES modelization consists in solving a part of this inertial range, in addition to the large scales of the flow. In particular, the Smagorinsky model is intrinsically discrete, linked to a space discretization that determines the large scales to solve. It adds an additional term to the averaged Navier-Stokes equations, called the eddy viscosity term, which is related to a grid Th that models the energy dissipated by the sub-grid (unresolved) scales. Then, the Smagorinsky turbulence model in its differential form is formulated as                    −1 Re∆w+(w·∇)w−∇·(νt(w)∇w)+∇p=f,in Ω, ∇·w=0,in Ω, w=g,on ΓD, −pn+1 Re +νt(w)∂w ∂n=0,on ΓN, (1.9) where Re is the Reynolds number defined in (A.22) , w is the velocity field, p is the pressure per mass density, νt(w) = ∑ K∈Th (CShK)2|∇w|K| is the eddy viscosity, where |· | denotes the Frobenius norm in Rd×d defined in (A.1) , the constant CS is the Smagorinsky constant which is estimated to be CS≈0.18 (cf. [ 31 , 26 , 44 ]) and hK is the diameter of the element K∈Th . In Chapter 6 in [ 25 ], V. John study the existence and uniqueness of solutions of the Smagorinsky model (1.9) following the results presented by Ladyzhenskaya in [ 30 ]. 20 1.4. ROM application to steady Smagorinsky turbulence model 1.4.1 Reduced Basis steady Smagorinsky model Let us state the basics of the RB approximation to the steady Smagorinsky model, considering the Reynolds number as the parameter µ in a suitable parameter set D . In order to apply the RB Method, we need a model with homogeneous Dirichlet BCs. To solve this, we decompose the velocity field into two fields, u the solution of a problem with homogeneous BCs and the lift uD , which should achieve the non-homogeneous Dirichlet BC uD=g on ΓD . Without lost of generality, we suppose that the lift is divergence free since this simplifies the problem. Finite Element approximation of Steady Smagorinsky model As in Section 1.1, we shall start from a variational formulation of our problem in a suitable Hilbert space and then build a discretization. Let us consider the spaces Y={v∈H1(Ω):v|ΓD=0},Q={q∈L2(Ω):ZΩ q=0} for velocity and pressure, and denote X=Y×Q . Let Yh and Qh be Finite Element spaces constructed on the grid Th to approximate Y and Q , respectively and Xh=Yh×Qh . Then the variational discretization of problem (1.9) is stated as follows:                  Find (uh,ph) = (uh(µ),ph(µ)) ∈Xhsuch that a(uh,vh;µ)+at(uh+uD;uh+uD,vh;µ) +b(vh,ph;µ)+c(uh,uh,vh;µ)+c(uD,uh,vh;µ) +c(uh,uD,vh;µ) = F(vh;µ),∀vh∈Yh, b(uh,qh;µ) = 0,∀qh∈Qh. (1.10) where the bilinear forms a(·,·;µ) and b(·,·;µ) are defined by a(u,v;µ) = 1 µZΩ ∇u:∇vdΩ,b(v,q;µ) = −ZΩ(∇·v)q dΩ; the trilinear form, c(·,·,·;µ) and the eddy diffusion term, at(·;·,·;µ) , are given by c(u,v,w;µ) = ZΩ(u·∇)v·wdΩ, at(u;v,w;µ) = ZΩ νt(u)∇v:∇wdΩ. Also, the linear form F(·;µ) is defined by F(v;µ) = ⟨f,v⟩Ω−a(uD,v;µ)−c(uD,uD,v;µ). 1. Basics of Reduced Basis modelling 21 Lastly, let us define the norms for the velocity and pressure spaces, considering the natural norm for each case. For the velocity space Y and, by extension, for Yh , we define the inner product as follows (u,v)T=ZΩ1 ˜ µ+˜ νt∇u:∇vdΩ,∀u,v∈Y(1.11) where ˜ νt=νt(uh(˜ µ) + uD) and ˜ µ=argmax µ∈D(∑ K∈Th (CShK)2min x∈K|∇uh(µ)+∇uD|(x)). This inner product induces a “turbulence” norm ∥·∥T= (·,·)1/2 T , which is equivalent to the H1 0(Ω) -norm. We use this norm since it takes into account the turbulence viscosity effects that governs the flow. For the pressure space Q and, by extension, for Qh , we consider the usual L2 -norm. Finally, for the X and Xh space, we consider the norm ∥V∥X=q∥v∥2 T+∥q∥2 L2(Ω),∀V= (v,q)∈X.(1.12) Following the notation used in (1.5) , we can state the discrete problem (Find Uh(µ) = (uh(µ),ph(µ)) ∈Xhsuch that A(Uh(µ),Vh;µ) = F(Vh;µ),∀Vh= (vh,qh)∈Xh (1.13) where A(Uh,Vh;µ) = 1 µA0(Uh,Vh)+A1(Uh,Vh)+A2(Uh,Vh)+A3(Uh,Vh)(1.14) and A0(U,V) = ZΩ ∇u:∇vdΩ, A1(U,V) = ZΩ(∇·u)q−(∇·v)p dΩ+ZΩ(uD·∇)u·v+(u·∇)uD·vdΩ, A2(U,V) = ZΩ(u·∇)u·vdΩ, A3(U,V) = ZΩ νt(u+uD)∇(u+uD):∇vdΩ. Note that now A(·,·;µ) is not bilinear, due to the convection and the eddy diffusion terms. Problem (1.13) is well posed, we will see it by means of the BRR theory of approximation of branches of regular solutions of non-linear problems in Section 1.4.2. 22 1.4. ROM application to steady Smagorinsky turbulence model RB Smagorinsky model Now, we are able to define the RB Problem for the Smagorinsky model. We define XN as the reduced space, then the RB problem is stated as (find UN(µ) = (uN(µ),pN(µ)) ∈XNsuch that A(UN(µ),VN;µ) = F(VN;µ),∀VN= (vN,qN)∈XN (1.15) The reduced space is XN=YN×QN , where YN and QN are the reduced spaces for velocity and pressure, respectively. The inf-sup condition is necessary to the well-posedness of problem (1.15) . Since YN is formed by solutions of (1.13) and these are nearly divergence-free, the pair velocity-pressure (YN,QN) is not inf-sup stable. In [ 11 ], it is proposed the use of an inner pressure supremizer operator Tµ N:Qh−→Yh defined as follows (Tµ Nph,vh)T=b(vh,ph;µ),∀vh∈Yh,(1.16) to enrich the velocity space YN in a Greedy Algorithm (cf. [ 47 ]). This ensures that YN , enriched with the supremizer, besides QN are inf-sup stable. The Greedy algorithm to construct the RB approximation of Smagorinsky model is described in Algorithm 6. Algorithm 6 (Weak) Greedy for steady Smagorinsky model Set εtol >0, Nmax ∈Nand S1={µ1}; Compute Uh(µ1) = (uh(µ1),ph(µ1)); Q1=span{ψ1:=ph(µ1)}; Y1=span{ξ ξ ξ1:=uh(µ1),ξ ξ ξ2:=Tµ Nψ1}; for N=2 : Nmax do µN=argmaxµ∈Dtrain ∆N−1(µ); εN−1=∆N−1(µ µ µN); if εN−1≤εtol then Nmax =N−1; end if Compute Uh(µN) = (uh(µN),ph(µN)); SN=SN−1∪{µN}; QN=QN−1∪span{ψN:=ph(µN)}; YN=YN−1∪span{ξ ξ ξ2N−1:=uh(µN),ξ ξ ξ2N:=Tµ NψN}; end for Remark 1.2. The use of an inner pressure supremizer Tµ N affects the dimension of the reduced spaces YN and QN . In our case, using the Greedy Algorithm, the dimension of the velocity space is twice that of the pressure space. 1. Basics of Reduced Basis modelling 23 Finally, we are able to define the reduced parameter independent matrices and tensors that will be stored in the offline phase. Then, we define AN,DN∈M2N×2N(R) , BN∈MN×2N(R) , CN∈M2N×2N×2N(R) defined by (AN)i j =RΩ∇ξ ξ ξj:∇ξ ξ ξidΩ,i,j=1,...,2N, (BN)i j =−RΩ(∇·ξ ξ ξj)ψidΩ,i=1,...,N,j=1,...,2N, (CN)i js =RΩ(ξ ξ ξs·∇)ξ ξ ξj·ξ ξ ξidΩ,i,j,s=1,...,2N, (DN)i j =RΩ(uD·∇)ξ ξ ξj·ξ ξ ξi+(ξ ξ ξj·∇)uD·ξ ξ ξidΩ,i,j=1,...,2N. where YN=span{ξ ξ ξi}2N i=1 and QN=span{ψi}N i=1 . Since the Smagorinsky term is nonlinear with respect to the parameter due to the eddy diffusion term, we use the Empirical Interpolation Method defined in Section 1.3 to approximate it. With this method, we get the RB space GM=span{qm}M m=1 and we obtain the following linear decomposition, νt(u)≈ M ∑ m=1 σj(µ)qm(x).(1.17) With the computation of the coefficients {σj(µ)}M j=1 for a given µ∈D , we are able to approximate the eddy diffusion term during the online phase. This yields the approximation at(u;u,v;µ µ µ)≈ M ∑ m=1 σj(µ)ZΩ qm∇u:∇vdΩ. Finally, we are able to define the matrix Em N∈M2N×2N(R) for m=1,...,M associated to the eddy diffusion term by (Em N)i j =ZΩ qm∇ξ ξ ξj:∇ξ ξ ξidΩ,i,j=1,...,2N. 1.4.2 A posteriori error bound estimator In this Section, we summarize Sections 4 and 5 in [ 11 ], where the a posteriori error bound estimator is derived. It is based on the Brezzi-Rappaz-Raviart (BRR) theory mentioned above (cf. [ 8 ]). Thanks to the BRR theory, it is proved that in a neighbourhood of UN(µ) solution of (1.15) , there exists a unique Uh(µ) solution of (1.13) , with ∥Uh(µ)−UN(µ)∥X≤∆N(µ),∀µ∈D 24 1.5. Application to 3D lid-driven cavity turbulent flow where ∆N(µ) = βN(µ) 2ρT(µ)h1−p1−τN(µ)i,τN(µ) = 4εN(µ)ρT(µ) βN(µ)2.(1.18) Here, •βN(µ)is the Stability Factor βN(µ)≡inf Zh∈Xh sup Vh∈Xh ∂1A(UN(µ µ µ),Vh;µ)(Zh) ∥Zh∥X∥Vh∥X ,(1.19) where the derivative operator of A(·,·;µ)follows the Definition A.1. •ρT is the Lipschitz continuity of ∂1A(UN(µ µ µ),Vh;µ)(Zh) for all Vh,Zh∈Xh defined by ρT=2C4;T+4C3;2CSh2−d/3 max where C4;T>0 is the Sobolev embedding constant verifying that ∥·∥L4(Ω)≤ C4;T∥·∥T,hmax =maxK∈ThhKand C3;2 >0 is defined in (A.11). •εN(µ)is the residual defined by εN(µ) = ∥R(UN(µ);µ)∥X′is the residual where ⟨R(Zh;µ),Vh⟩=A(Zh,Vh;µ)−F(Vh;µ),∀Zh,Vh∈Xh. Note that (1.18) only applies when the dual norm of the residual εN(µ) is small enough to have τ(µ)N<1 for all µ∈D . 1.5 Application to 3D lid-driven cavity turbulent flow In this Section, we extend Section 7.2 in [ 11 ] to the 3D case, actually to 3D lid-driven cavity flow, using the same software, FreeFem++ v. 4.8 (cf. [ 23 ]). To solve the problem presented, we use a FE approximation with the Taylor-Hood finite element, i.e., we consider P2−P1 for velocity-pressure discretization that are inf-sup stable. Since the pressure is in L2(Ω) with a mean of zero, we add a stabilization term for the pressure discretization in the variational formulation, this is, a L2 penalization. We consider a lid-driven cavity flow problem in the unit cube as we see in Figure 1.1a. We define a tetrahedral mesh in the cube, with 13 subdivisions for each edge, obtaining a mesh with 13182 tetrahedral and 2744 vertices (see Figure 1.1b). Furthermore, the grid element diameter is hK=0.1332 for each K∈Th . The parameter set is an interval such that D= [Remin,Remax] with Remin =1000 and Remax =5000. In this range, the flow is steady and laminar, therefore, this test is purely 1. Basics of Reduced Basis modelling 25 yx z ΓS ΓW (a) 3D domain. (b) Mesh. Figure 1.1: Unit cube for the problem. academic to verify the estimate (1.18) . Recall that the eddy viscosity term adds diffusion to our problem and thanks to this, we are able to achieve a steady situation. As we consider a lid-driven cavity problem, we only need to set Dirichlet BCs. We identify two parts on the boundary, the top ΓS and the walls ΓW as it is shown in Figure 1.1a. We impose that w|ΓW= 0and w|ΓS= (g(x,y),0) with g(x,y) = 16x(x−1)y(y−1) for all (x,y)∈ΓS . As we already mentioned, we need to split the velocity field w into u and a lift uD . For this case, we choose uD to be the solution of Stokes problem. Since the problem (1.10) is nonlinear, we use a semi-implicit linearized evolution approach,                              Assuming that u0 h=0, then for k=1,...,NFE , find (uk h,pk h) = (uk h(µ µ µ),pk h(µ µ µ)) ∈Xhsuch that 1 ∆t(uk h,vh)Ω+a(uk h,vh;µ)+at(uk−1 h+uD;uk h+uD,vh;µ) +b(vh,pk h;µ)+c(uk−1 h,uk h,vh;µ)+c(uD,uk h,vh;µ) +c(uk h,uD,vh;µ) = F(vh;µ)+ 1 ∆t(uk−1 h,vh)Ω,∀vh∈Yh, b(uk h,qh;µ) = 0,∀qh∈Qh. (1.20) We stop the process when a nearly steady solution has been reached, this is, when for k=1,...,NFE we have that ∥(uk h−uk−1 h)/∆t∥L2(Ω)<εFE for εFE =10−11 and ∆t=1 . The selection to εFE affects directly to the residual εN(µ) for all µ∈D in 32 Part I Multiparameter Reduced Basis steady Smagorinsky model 2 Reduced Basis steady Smagorinsky model with variable geometry In this chapter, we build a RB model for the air flow within a cloister, as a problem of practical interest of a turbulent flow with multi-parametric geometrical dependence. We consider a steady turbulent flow governed by the 2D Smagorinsky model. The geometrical parameters are the width and the height of the corridors around the cloister. In Section 2.1 we define the problem establishing the domain, the geometrical parameters, and the procedure to follow for applying the RB Method. In Section 2.2 we define the RB problem explaining how the RB matrices are built. Finally, in Section 2.3 we deduce an a posteriori error estimator necessary to build the RB space in the Greedy Algorithm. At the end of the chapter, in Section 2.4, we show the numerical results for a specific case, giving keys for coding and validating the estimator developed in Section 2.3. 2.1 Problem statement We consider a 2D modeling of the air flow within a cloister. Although this will avoid 3D effects in the aerothermal flow, this still is an adequate approach to determine its optimal dimensions regarding thermal comfort. Indeed, thermal inversions and strong 3D effects will only be produced by strong winds that are assumed to occur only occasionally. We thus consider thermally stratified flows through a 2D modeling. 35 36 2.1. Problem statement We consider the 2D steady Smagorinsky model that allows us to solve the largest scales considering a coarse mesh. This model was introduced in Section 1.4. However, we will derive the dimensionless model as we depart from a real situation. The flow domain is described in Figure 2.1 which can be considered as a vertical section of the real domain (a cloister), where ω and σ are the width and the height of the corridor around the cloister, respectively, and W and H are the total width and height. The parameters for this case are µ µ µ= (ω,σ) while W and H are fixed. We denote by Ω(µ µ µ)∈R2 the domain for each µ µ µ∈D with D the parameter set. We take D= [ωmin,ωmax]×[σmin,σmax] where ωmin , ωmax , σmin and σmax are the extreme values of corridors width and height that make sense for architectural design. In the cases σ=0 and σ=H , there would be no cloister, it would be a question of studying a courtyard, then [σmin,σmax]⊂(0,H) . For ω , we consider only up to W/2 due to symmetry and we exclude the extremes for the same reason as for the height, then [ωmin,ωmax]⊂(0,W/2) . (0,0) (0,σ) (ω,σ) (ω, H)(W−ω, H) (W−ω,σ) (W, σ) (W, 0) u|ΓF= (0,0) u|ΓS= (g, 0) u|ΓW= (0,0) (ω,0) Figure 2.1: Domain and diagram of the problem. Let Γ(µ µ µ) = ∂Ω(µ µ µ) be the domain boundary that splits into three parts: • Open sky boundary: ΓS(µ µ µ), • Floor boundary: ΓF(µ µ µ), • Lateral walls boundary: ΓW(µ µ µ). 2. RB steady Smagorinsky model with variable geometry 37 The flow within the cloister is induced through boundary ΓS . This setup leads us to a lid-driven cavity turbulent flow problem. We impose homogeneous Dirichlet boundary conditions (BCs) on the velocity on the floor ΓF(µ µ µ) and the walls ΓW(µ µ µ) . On the top we assume an horizontal wind, given by the function g(x1)≥0 with x1∈ΓS . We set {Th}h>0 as a family of triangulations of Ω . In this way, the problem to solve is: Find w:Ω(µ µ µ)−→ R2 and p:Ω(µ µ µ)−→ R such that            −ν∆w+(w·∇)w−∇·(νt(w;µ µ µ)∇w)+∇p=0,in Ω(µ µ µ), ∇·w=0,in Ω(µ µ µ), w= (g(·;µ µ µ),0),on ΓS(µ µ µ), w=0,on ΓW(µ µ µ)∪ΓF(µ µ µ), (2.1) where w is the velocity field, p is the pressure per mass density, ν is the kinematic viscosity for the air and the turbulent viscosity is given by νt(w;µ µ µ) = ∑ K∈Th (CShK(µ µ µ))2|∇w(µ µ µ)|K|, where we recall that hK is the diameter of the element K∈Th and CS denotes the Smagorinsky constant. As it is usual, we consider the dimensionless model. In this way, we use a single model for flows within similar domains, and deal with computer quantities close to one to avoid round-off errors. We set U0 and L0 as the characteristics velocity and length, respectively, then w=U0w∗ , p=ρU2 0p∗ , µ µ µ=L0µ µ µ∗ , where w∗ and p∗ represent the dimensionless velocity and pressure, µ µ µ∗ the dimensionless length parameter, ρ is the mass density of the air (constant in this case), g=U0g∗ , hK=L0h∗ K . Henceforth, (w∗,p∗) satisfy the problem                −1 Re∆∗w∗+(w∗·∇∗)w∗−∇∗·(ν∗ t(w∗;µ µ µ∗)∇∗w∗)+∇∗p∗=0,in Ω(µ µ µ∗), ∇∗·w∗=0,in Ω(µ µ µ∗), w∗= (g∗(·;µ µ µ∗),0),on ΓS(µ µ µ∗), w∗=0,on ΓW,F(µ µ µ∗), where Re =U0L0/ν is the Reynolds number defined in (A.22) , fixed for this problem, and ν∗ t(w∗;µ µ µ∗) = ∑ K∈Th (CSh∗ K(µ µ µ∗))2|∇∗w∗(µ µ µ∗)|.(2.2) The operator ∇∗ and ∆∗ are defined respect to the dimensionless variable x∗=x/L0 . From now on, we avoid the star notation for simplicity. 38 2.1. Problem statement In order to apply the RB Method, we need a model with homogeneous Dirichlet BCs. To solve this, we decompose the velocity field as the sum of two fields, u the solution of a problem with homogeneous Dirichlet BCs and uD a function defined in all Ω(µ µ µ) sufficiently regular and that satisfies non-homogeneous Dirichlet BCs, this is, •uD= (g(·;µ µ µ),0)on ΓS(µ µ µ), •uD=0on ΓW(µ µ µ)∪ΓF(µ µ µ). This function is known as the boundary lifting function or simply a lift of the conditions. Then, the problem leads For a known uD , find (u,p) such that                  −1 Re∆u+(u·∇)u+(uD·∇)u+(u·∇)uD+∇p −∇·(νt(u+uD;µ µ µ)∇(u+uD)) = 1 Re∆uD−(uD·∇)uD,in Ω(µ µ µ), ∇·u=−∇·uD,in Ω(µ µ µ), u=0,on Γ(µ µ µ). (2.3) The total solution w=u+uD is completely independent of the lift choice. This uD can be determined by an analytical expression or by the solution of a simpler problem which preserves the nonhomogeneous BCs, for example, the Stokes problem. We are interested in looking at weak solutions to establish an approximation, and therefore we start by defining the variational formulation for this problem. Let us consider the spaces Y(µ µ µ) = H1 0(Ω(µ µ µ)) , Q(µ µ µ) = L2 0(Ω(µ µ µ)) , and X(µ µ µ) = Y(µ µ µ)×Q(µ µ µ) . We can define the weak formulation of the problem (2.3) , doing the standard procedure:                        Find (u,p) = (u(µ µ µ),p(µ µ µ)) ∈X(µ µ µ)such that a(u,v;µ µ µ)+b(v,p;µ µ µ)+c(u,u,v;µ µ µ) +c(uD,u,v;µ µ µ)+c(u,uD,v;µ µ µ) +at(u+uD;u+uD,v;µ µ µ) =−a(uD,v;µ µ µ)−c(uD,uD,v;µ µ µ),∀v∈Y(µ µ µ), b(u,q;µ µ µ) = −b(uD,q;µ µ µ),∀q∈Q(µ µ µ), (2.4) where the bilinear forms a(·,·;µ µ µ) and b(·,·;µ µ µ) are given by a(u,v;µ µ µ) = 1 Re ZΩ(µ µ µ) ∇u:∇vdΩ,b(v,q;µ µ µ) = −ZΩ(µ µ µ)(∇·v)q dΩ; 2. RB steady Smagorinsky model with variable geometry 39 the trilinear form, c(·,·,·;µ µ µ) and the eddy viscosity term, at(·;·,·;µ µ µ) , are given by c(u,v,z;µ µ µ) = ZΩ(µ µ µ)(u·∇)v·zdΩ, at(u;v,z;µ µ µ) = ZΩ(µ µ µ) νt(u;µ µ µ)∇v:∇zdΩ. Ω1 Ω4 Ω2Ω3 Figure 2.2: Partitioned Geometry. These forms depend on the geometrical parameters that determine the domain. To apply the RB Method, we need to express these forms as linear expressions with respect to the parameters. With this purpose, we consider a reference domain, on which we build the variational formulation through a change of variables. We set µ µ µr= (ωr,σr)∈D a couple of reference parameters that define Ωr=Ω(µ µ µr) as the reference domain. Then, we consider a smooth mapping Φ(·;µ µ µ) that transforms the reference domain into the original domain for any value of the parameter µ µ µ . We have defined Φas a piecewise linear mapping and because of that we divide the domain in four regions Ωl for l=1,2,3,4 as is shown in Figure 2.2 and we define a linear transformation for each one of them, this is, Φl for l=1,2,3,4 such that Φ(x;µ µ µ) = 40 2.1. Problem statement Φl(x;µ µ µ) if x∈Ωl , as follows Φ1(x;µ µ µ) =    ω ωr 0 0σ σr    x1 x2!, Φ2(x;µ µ µ) =     W−2ω W−2ωr 0 0σ σr     x1 x2!+  ω−W−2ω W−2ωr 0 , Φ3(x;µ µ µ) =    ω ωr 0 0σ σr    x1 x2!+  Wωr−ω ωr 0 , Φ4(x;µ µ µ) =     W−2ω W−2ωr 0 0H−σ H−σr     x1 x2!+    ω−ωr W−2ω W−2ωr σ−σr H−σ H−σr     . (2.5) Observe that the global mapping is Lipschitz-continuous as Φ1|Γ12 =Φ2|Γ12,Φ2|Γ23 =Φ3|Γ23 and Φ2|Γ24 =Φ4|Γ24, where Γi j =∂Ωi∩Ωj ; and Φl∈C∞(Ωl) for l=1,2,3,4 . The Jacobian determinant for each transformation is |JΦ1|=|JΦ3|=ωσ ωrσr ,|JΦ2|=(W−2ω)σ (W−2ωr)σr ,|JΦ4|=(W−2ω)(H−σ) (W−2ωr)(H−σr),(2.6) Associated with the Jacobian matrix for each transformation, we define three diagonal matrices ηl=ηl(µ µ µ),λl=λl(µ µ µ),ϕl=ϕl(µ µ µ)∈M2×2(R) for l=1,2,3,4 by ηl=J−T Φl|JΦl|,λl=J−1 ΦlJ−T Φl|JΦl|,ϕl=J−1 ΦlJ−T Φl, whose diagonal elements are: η1 11 =η2 11 =η3 11 =σ σr ,η4 11 =H−σ H−σr ,η1 22 =η3 22 =ω ωr ,η2 22 =η4 22 =W−2ω W−2ωr ,(2.7) λ1 11 =λ3 11 =ωrσ ωσr ,λ2 11 =(W−2ωr)σ (W−2ω)σr ,λ4 11 =(W−2ωr)(H−σ) (W−2ω)(H−σr),(2.8) λ1 22 =λ3 22 =ωσr ωrσ,λ2 22 =(W−2ω)σr (W−2ωr)σ,λ4 22 =(W−2ω)(H−σr) (W−2ωr)(H−σ),(2.9) 2. RB steady Smagorinsky model with variable geometry 41 and ϕ1 11 =ϕ3 11 =ω2 r ω2,ϕ2 11 =ϕ4 11 =(W−2ωr)2 (W−2ω)2,(2.10) ϕ1 22 =ϕ2 22 =ϕ3 22 =σ2 r σ2,ϕ4 22 =(H−σr)2 (H−σ)2.(2.11) These matrices are well defined since µ µ µ∈D⊂(0,W/2)×(0,H) . Now, we can state the weak formulation (2.4) on Ωr . Let us consider the spaces Yr=Y(µ µ µr) , Qr=Q(µ µ µr) , and Xr=Yr×Qr . Let us define the norms related to the spaces Yr and Qr . For the velocity space Yr , we consider the turbulence norm which is the natural norm associated to the problem, defined as (v,z)T=ZΩr1 Re +˜ νt∇v:∇zdΩ,∀v,z∈Yr(2.12) where ˜ νt=νt(u+uD;˜ µ µ µ) and ˜ µ µ µ=arg min µ µ µ∈D∑ K∈Th (CShK(µ µ µ))2min x∈K|∇(u(µ µ µ)+ uD)|(x). The inner product (·,·)T is parameter independent, it is well defined and it induces a norm ∥· ∥T= (·,·)1/2 T in Yr equivalent to the H1 0 -norm. For the pressure space Qr , we will use the usual L2 -norm. Finally, we define the X -norm such as ∥V∥X=q∥v∥2 T+∥q∥2 L2(Ω),∀V= (v,q)∈Xr.(2.13) Then, the weak formulation reads:                                  Find (u,p) = (u(µ µ µ),p(µ µ µ)) ∈Xrsuch that 4 ∑ l=1 2 ∑ k=1 λl kk(akl(u,v)+at,kl(u+uD;u+uD,v;µ µ µ)) +ηl kk(bkl(v,p)+ckl(u,u,v)+ckl(uD,u,v)+ckl(u,uD,v)) =− 4 ∑ l=1 2 ∑ k=1 λl kkakl(uD,v)+ηl kkckl(uD,uD,v),∀v∈Yr, 4 ∑ l=1 2 ∑ k=1 ηl kkbkl(u,q) = − 4 ∑ l=1 2 ∑ k=1 ηl kkbkl(uD,q),∀q∈Qr. (2.14) Notice that the indices kand lare associated to the derivative direction and the different regions in Ωr , respectively. Then, for k=1,2 and l=1,2,3,4 the forms associated to the problem (2.14) are defined by akl(u,v) = ZΩl 1 Re(∂ku1∂kv1+∂ku2∂kv2)dΩ,(2.15) 48 2.3. A posteriori error estimator Then, we use (2.29) and (2.30) and the definition of the X-norm (2.13), |∂1A1(Uh,Vh;µ µ µ)(Zh)|≤η∗C1,2;T(∥vh∥T∥rh∥L2(Ω)+∥zh∥T∥qh∥L2(Ω)) +η∗C4;T(C1,2;T∥uD∥L4(Ω)+C4;T∥∇uD∥L2(Ω))∥vh∥T∥zh∥T ≤η∗(C1,2;T+C4;TC)∥Vh∥X∥Zh∥X, where C= (C1,2;T∥uD∥L4(Ω)+C4;T∥∇uD∥L2(Ω)). • The following term, |∂1A2(Uh,Vh;µ µ µ)(Zh)|= 4 ∑ l=1 2 ∑ k=1 ηl kk(ckl(uh,zh,vh)+ckl(zh,uhvh)) ≤η∗ZΩ(uh·∇)zh·vh+(zh·∇)uh·vhdΩ ≤η∗∥uh∥L4(Ω)∥∇zh∥L2(Ω)∥vh∥L4(Ω) +∥zh∥L4(Ω)∥∇uh∥L2(Ω)∥vh∥L4(Ω)) ≤η∗C1,2;TC4;T∥uh∥L4(Ω)+C2 4;T∥∇uh∥L2(Ω)∥vh∥T∥zh∥T ≤η∗C1,2;TC4;T∥uh∥L4(Ω)+C2 4;T∥∇uh∥L2(Ω)∥Vh∥X∥Zh∥X. • Finally, |∂1A3(Uh,Vh,µ µ µ)(Zh)|≤ 4 ∑ l=1 2 ∑ k=1 λl kkat,kl(wh;zh,vh;µ µ µ) + 4 ∑ l=1∑ K∈Th,lZK (CShK)2∑2 k=1λl kk(∇kwh·∇kzh) q∑2 k=1ϕl kk|∇kwh| 2 ∑ k=1 ϕl kk(∇kwh·∇kvh)! ≤λ∗∑ K∈ThZK (CShK(µ µ µ))2|∇wh|∇zh:∇vh +λ∗ϕ∗ pˆ ϕ∑ K∈ThZK (CShK(µ µ µ))2∇wh:∇zh |∇wh|(∇wh:∇vh) ≤λ∗ 1+ϕ∗ pˆ ϕ!∑ K∈ThZK (CShK(µ µ µ))2|∇wh||∇zh||∇vh| 2. RB steady Smagorinsky model with variable geometry 49 ≤λ∗ 1+ϕ∗ pˆ ϕ!∑ K∈Th (CShK(µ µ µ))2∥∇wh∥L3(K)∥∇zh∥L3(K)∥∇vh∥L3(K) ≤λ∗ 1+ϕ∗ pˆ ϕ!C2 3;2 ∑ K∈Th C2 ShK(µ µ µ)2−d/3∥∇wh∥L3(K)∥∇zh∥L2(K)∥∇vh∥L2(K) ≤λ∗ 1+ϕ∗ pˆ ϕ!C2 S(h∗ max)2−d/3C2 3;2∥∇wh∥L3(Ω)∥∇zh∥L2(Ω)∥∇vh∥L2(Ω) ≤λ∗ 1+ϕ∗ pˆ ϕ!C2 S(h∗ max)2−d/3C2 1,2;TC2 3;2∥∇wh∥L3(Ω)∥zh∥T∥vh∥T ≤λ∗ 1+ϕ∗ pˆ ϕ!C2 S(h∗ max)2−d/3C2 1,2;TC2 3;2∥∇wh∥L3(Ω)∥Vh∥X∥Vh∥X. where we have used (2.31) and that wh=uh+uD. Consequently, the continuity constant is γh(µ µ µ) = λ∗C2 1,2;T Re +η∗(C1,2;T+C4;TC) +η∗C1,2;TC4;T∥uh∥L4(Ω)+C2 4;T∥∇uh∥L2(Ω) +λ∗ 1+ϕ∗ pˆ ϕ!C2 S(h∗ max)2−d/3C2 1,2;TC2 3;2∥∇uh+∇uD∥L3(Ω). where the dependency on µ µ µis produced by uh. Now, we study the coercivity of the derivative operator. Proposition 2.2. For any µ µ µ∈D and Uh=Uh(µ µ µ)∈Xh , let us suppose that ∥∇uD∥L2(Ω)<1/C and ∥∇uh∥L2(Ω)≤1 C− ∥∇uD∥L2(Ω), where the positive constant C is defined by C= ˆ λmin{1,pˆ ϕ} 2η∗C2 1,2;TC2 4;1,2 with C1,2;Tand C4;1,2defined in (2.29)-(2.30). Then it holds, ∂1A(Uh,Zh;µ µ µ)(Zh)≥˜ βh∥zh∥2 T,∀Zh∈Xh, 50 2.3. A posteriori error estimator where ˜ βh=ˆ λmin{1,pˆ ϕ}−2η∗C2 1,2;TC2 4;1,2(∥∇uh∥L2(Ω)+∥∇uD∥L2(Ω))>0. Proof. By definition, ∂1A(Uh,Zh;µ µ µ)(Zh) = 3 ∑ i=0 ∂1Ai(Uh,Zh;µ µ µ)(Zh) and we study each term separately. We recall that forms Ai for i=0,1,2,3 are defined in (2.21). Along the proof, we will use that wh=uh+uDfor simplicity. • The derivative of the diffusion term is bounded as ∂1A0(Uh,Zh;µ µ µ)(Zh) = 4 ∑ l=1 2 ∑ k=1 λl kkakl(zh,zh)≥ˆ λ1 Re∥∇zh∥2 L2(Ω). The eddy viscosity term is treated by ∂1A3(Uh,Zh;µ µ µ)(Zh) = 4 ∑ l=1 2 ∑ k=1 λl kkat,kl(wh;zh,zh;µ µ µ) + 4 ∑ l=1∑ K∈Th,lZK (CShK)2 2 ∑ k=1 λl kk∇kwh·∇kzh s2 ∑ k=1 ϕl kk|∇kwh|2 2 ∑ k=1 ϕl kk∇kwh·∇kzh! ≥ˆ λpˆ ϕ∑ K∈ThZK (CShK)2|∇wh||∇zh|2+ ˆ λˆ ϕ √ϕ∗∑ K∈ThZK (CShK)2|∇wh:∇zh|2 |∇wh|. Using that this last addend is positive and adding these two terms, we obtain ∂1A0(Uh,Zh;µ µ µ)(Zh)+∂1A3(Uh,Zh;µ µ µ)(Zh) ≥ˆ λ1 Re ZΩ|∇zh|2dΩ+ˆ λpˆ ϕZΩ(CShK)2|∇wh||∇zh|2dΩ + ˆ λˆ ϕ √ϕ∗∑ K∈ThZK (CShK)2|∇wh:∇zh|2 |∇wh| ≥ˆ λmin{1,pˆ ϕ}∥zh∥2 T. (2.32) 2. RB steady Smagorinsky model with variable geometry 51 • On the other hand, we bound from above the rest of the terms, such as |∂1A1(Uh,Zh;µ µ µ)(Zh)+∂1A2(Uh,Zh;µ µ µ)(Zh)| = 4 ∑ l=1 2 ∑ k=1 ηl kk(ckl(wh,zh,zh)+ckl(zh,wh,zh)) ≤η∗ZΩ(wh·∇)zh·zh+(zh·∇)wh·zhdΩ ≤η∗ZΩ|wh||∇zh||zh|+|zh|2|∇wh|dΩ ≤η∗∥wh∥L4(Ω)∥∇zh∥L2(Ω)∥zh∥L4(Ω)+∥zh∥2 L4(Ω)∥∇wh∥L2(Ω) ≤2η∗C2 4;1,2∥∇wh∥L2(Ω)∥∇zh∥2 L2(Ω) ≤2η∗C2 1,2;TC2 4;1,2∥∇wh∥L2(Ω)∥zh∥2 T where C4;1,2 is defined in (2.30) and thanks to the Triangle Inequality, we obtain that |∂1A1(Uh,Zh;µ µ µ)(Zh)+∂1A2(Uh,Zh;µ µ µ)(Zh)| ≤2η∗C2 1,2;TC2 4;1,2(∥∇uh∥L2(Ω)+∥∇uD∥L2(Ω))∥zh∥2 T. Then, ∂1A1(Uh,Zh;µ µ µ)(Zh)+∂1A2(Uh,Zh;µ µ µ)(Zh) ≥−2η∗C2 1,2;TC2 4;1,2(∥∇uh∥L2(Ω)+∥∇uD∥L2(Ω))∥zh∥2 T.(2.33) Finally, regrouping (2.32) and (2.33), we obtain ∂1A(Uh,Zh;µ µ µ)(Zh)≥˜ βh∥zh∥2 T, where ˜ βh=ˆ λmin{1,pˆ ϕ}−2η∗C2 1,2;TC2 4;1,2(∥∇uh∥L2(Ω)+∥∇uD∥L2(Ω)) Since ˜ βhshould be positive, we define C>0 C= ˆ λmin{1,pˆ ϕ} 2η∗C2 1,2;TC2 4;1,2 then, if ∥∇uD∥L2(Ω)<1/C and ∥∇uh∥L2(Ω)≤1/C−∥∇uD∥L2(Ω) , we are able to guarantee that the constant is positive and therefore, the operator is coercive. 52 2.3. A posteriori error estimator Remark 2.3. If we choose Yh and Qh in order to verify the inf-sup condition in Theorem A.8 (for example, considering Taylor-Hood finite element spaces), we are able to ensure that the operator ∂1A verifies (2.28) . Moreover, if the pair of spaces (YN,QN) , are inf-sup stables and the data are small enough, the problem (2.26) is also well-posed. According to BRR theory, we can conclude that problem (2.22) admits a solution Uh∈Xh for small data which is locally unique. 2.3.2 The building of an a posteriori error estimator As a previous result, we need to prove that the tangent operator ∂1A is Lipschitz. Lemma 2.1. For all µ µ µ∈Dand U1 h,U2 h,Zh,Vh∈Xh, it holds |∂1A(U1 h,Vh;µ µ µ)(Zh)−∂1A(U2 h,Vh;µ µ µ)(Zh)|≤ρT∥U1 h−U2 h∥X∥Zh∥X∥Vh∥X(2.34) for a positive constant ρT=2η∗C1,2;TC2 4;T+ λ∗pϕ∗+3λ∗ϕ∗ pˆ ϕ!C2 S(h∗ max)2−d/2C3 1,2;TC3 3;2, where constant C1,2;T, C4;Tand C3;2 are defined in (2.29)-(2.31). Proof. By definition ∂1A(U1 h,Vh;µ µ µ)(Zh)−∂1A(U2 h,Vh;µ µ µ)(Zh) = ∂1A2(U1 h,Vh;µ µ µ)(Zh) −∂1A2(U2 h,Vh;µ µ µ)(Zh)+∂1A3(U1 h,Vh;µ µ µ)(Zh)−∂1A3(U2 h,Vh;µ µ µ)(Zh),(2.35) remaining the terms associated with the nonlinear part of the operator A . First, we study the first couple of terms above ∂1A2(U1 h,Vh;µ µ µ)(Zh)−∂1A2(U2 h,Vh;µ µ µ)(Zh) = 4 ∑ l=1 2 ∑ k=1 ηl kk(ckl(u1 h,zh,vh)−ckl(u2 h,zh,vh)+ckl(zh,u1 h,vh)−ckl(zh,u2 h,vh)) = 4 ∑ l=1 2 ∑ k=1 ηl kk(ckl(u1 h−u2 h,zh,vh)+ckl(zh,u1 h−u2 h,vh)) ≤η∗ZΩ(u1 h−u2 h)·∇zh·vhdΩ+ZΩ(zh·∇)(u1 h−u2 h)·vhdΩ. 2. RB steady Smagorinsky model with variable geometry 53 Then, ZΩ(u1 h−u2 h)·∇zh·vhdΩ≤ZΩ|u1 h−u2 h||∇zh||vh|dΩ ≤∥u1 h−u2 h∥L4(Ω)∥∇zh∥L2(Ω)∥vh∥L4(Ω)≤C1,2;TC2 4;T∥u1 h−u2 h∥T∥zh∥T∥vh∥T ≤C1,2;TC2 4;T∥U1 h−U2 h∥X∥Zh∥X∥Vh∥X, and analogously, ZΩ(zh·∇)(u1 h−u2 h)·vhdΩ≤∥zh∥L4(Ω)∥∇(u1 h−u2 h)∥L2(Ω)∥vh∥L4(Ω) ≤C1,2;TC2 4;T∥U1 h−U2 h∥X∥Zh∥X∥Vh∥X, and finally, |∂1A2(U1 h,Vh;µ µ µ)(Zh)−∂1A2(U2 h,Vh;µ µ µ)(Zh)| ≤2η∗C1,2;TC2 4;T∥U1 h−U2 h∥X∥Zh∥X∥Vh∥X,(2.36) where the constants C1,2;Tand C4;Twere defined in (2.29) and (2.30), respectively. Now we study the second couple of terms in (2.35), ∂1A3(U1 h,Vh;µ µ µ)(Zh)−∂1A3(U2 h,Vh;µ µ µ)(Zh) = 4 ∑ l=1 2 ∑ k=1 λl kk(at,kl(w1 h;zh,vh;µ µ µ)−at,kl(w2 h;zh,vh;µ µ µ)) + 4 ∑ l=1∑ K∈Th,lZK (CShK)2  ∑2 k=1λl kk(∇kw1 h·∇kzh) q∑2 k=1ϕl kk(∇kw1 h·∇kw1 h) 2 ∑ k=1 ϕl kk(∇kw1 h·∇kvh)! −∑2 k=1λl kk(∇kw2 h·∇kzh) q∑2 k=1ϕl kk(∇kw2 h·∇kw2 h) 2 ∑ k=1 ϕl kk(∇kw2 h·∇kvh)! . (2.37) where we recall that wi h=ui h+uD for i=1,2 and the form at,kl(·;·,·;µ µ µ) was defined in (2.18). It holds  4 ∑ l=1 2 ∑ k=1 λl kk(at,kl(w1 h;zh,vh;µ µ µ)−at,kl(w2 h;zh,vh;µ µ µ)) ≤λ∗ZΩ(νt(w1 h;µ µ µ)−νt(w2 h;µ µ µ))(∇zh:∇vh)dΩ ≤λ∗pϕ∗ ∑ K∈ThZK (CShK)2|∇w1 h|−|∇w2 h|(∇zh:∇vh) 54 2.3. A posteriori error estimator ≤λ∗pϕ∗∑ K∈ThZK (CShK)2|∇(u1 h−u2 h)||∇zh||∇vh| ≤λ∗pϕ∗∑ K∈Th (CShK)2∥∇(u1 h−u2 h)∥L3(K)∥∇zh∥L3(K)∥∇vh∥L3(K) ≤λ∗pϕ∗C3 3;2 ∑ K∈Th (CShK)2−d/2∥∇(u1 h−u2 h)∥L2(K)∥∇zh∥L2(K)∥∇vh∥L2(K) ≤λ∗pϕ∗C2 S(h∗ max)2−d/2C3 3;2∥∇(u1 h−u2 h)∥L2(Ω)∥∇zh∥L2(Ω)∥∇vh∥L2(Ω) ≤λ∗pϕ∗C2 S(h∗ max)2−d/2C3 1,2;TC3 3;2∥U1 h−U2 h∥X∥Zh∥X∥Vh∥X. where we have applied the Local Inverse Inequality (2.31). For the second term, using the same results and procedure that in the first term, we obtain:  4 ∑ l=1∑ K∈Th,lZK (CShK)2  ∑2 k=1λl kk(∇kw1 h·∇kzh) q∑2 k=1ϕl kk(∇kw1 h·∇kw1 h) 2 ∑ k=1 ϕl kk(∇kw1 h·∇kvh)! −∑2 k=1λl kk(∇kw2 h·∇kzh) q∑2 k=1ϕl kk(∇kw2 h·∇kw2 h) 2 ∑ k=1 ϕl kk(∇kw2 h·∇kvh)!  ≤λ∗ϕ∗ pˆ ϕ ∑ K∈ThZK (CShK)2∇w1 h:∇zh |∇w1 h|(∇w1 h:∇vh)−∇w2 h:∇zh |∇w2 h|(∇w2 h:∇vh) =λ∗ϕ∗ pˆ ϕ ∑ K∈ThZK (CShK)2∇w1 h:∇zh |∇w1 h|(∇(w1 h−w2 h):∇vh) +∇(w1 h−w2 h):∇zh |∇w2 h|(∇w2 h:∇vh) +∑ K∈ThZK (CShK)2(|∇w2 h|−|∇w1 h|) |∇w1 h||∇w2 h|(∇w1 h:∇zh)(∇w2 h:∇vh) ≤3λ∗ϕ∗ pˆ ϕ∑ K∈ThZK (CShK)2|∇(u1 h−u2 h)||∇zh||∇vh| ≤3λ∗ϕ∗ pˆ ϕ∑ K∈Th (CShK)2∥u1 h−u2 h∥L3(K)∥zh∥L3(K)∥vh∥L3(K) ≤3λ∗ϕ∗ pˆ ϕC2 S(h∗ max)2−d/2C3 1,2;TC3 3;2∥U1 h−U2 h∥X∥Zh∥X∥Vh∥X. Using these bounds in (2.37) and summing with (2.36) , we conclude that there exists ρT=2η∗C1,2;TC2 4;T+ λ∗pϕ∗+3λ∗ϕ∗ pˆ ϕ!C2 S(h∗ max)2−d/2C3 1,2;TC3 3;2, 2. RB steady Smagorinsky model with variable geometry 55 such that (2.34) is satisfied. Remark 2.4. Note that the Lipschitz constant ρT only depends on possible values of µ µ µ , this is, the selection of Dand the characteristic mesh parameters. It also depends on (h∗ max)2−d/2which tends to zero since 2−d/2>0for any d <4. As we explained, in order to apply the RB Method, we develop an a posteriori error estimator. For its construction, we first introduce some definitions. Definition 2.1. The Stability Factor βN(µ µ µ) and the continuity factor γN(µ µ µ) , associated with ∂1A(UN(µ µ µ),·;µ µ µ)are defined by βN(µ µ µ)≡inf Zh∈Xh sup Vh∈Xh ∂1A(UN(µ µ µ),Vh;µ µ µ)(Zh) ∥Zh∥X∥Vh∥X ,(2.38) γN(µ µ µ)≡sup Zh∈Xh sup Vh∈Xh ∂1A(UN(µ µ µ),Vh;µ µ µ)(Zh) ∥Zh∥X∥Vh∥X ,(2.39) for a given UN(µ µ µ)∈XN, solution of (2.26). Note that from Proposition 2.2, βN(µ µ µ)>0 whenever the data is small enough and the pair (Yh,Qh) are inf-sup stable. This condition implies that the tangent operator ∂1A(UN(µ µ µ),·;µ µ µ) is an isomorphism of Xh into its self for all µ µ µ∈D . Definition 2.2. For a given UN(µ µ µ)∈XN solution of (2.26) , the supremizer operator Tµ µ µ N:Xh−→Xhis defined by Tµ µ µ NZh=arg sup Vh∈Xh ∂1A(UN(µ µ µ),Vh;µ µ µ)(Zh) ∥Vh∥X .(2.40) This supremizer is unlike the inner-pressure supremizer presented in (2.23) . The pressure supremizer is related to the pressure and it is defined to enrich the velocity space, while this last supremizer is related to the tangent operator ∂1A . Remark 2.5. From definitions 2.1 and 2.2, we deduce the relation between the factors and the supremizer operator: βN(µ µ µ)≡inf Zh∈Xh ∥Tµ µ µ NZh∥X ∥Zh∥X ,γN(µ µ µ)≡sup Zh∈Xh ∥Tµ µ µ NZh∥X ∥Zh∥X . In particular, we have βN(µ µ µ)≤∂1A(UN(µ µ µ),Tµ µ µ NZh;µ µ µ)(Zh) ∥Tµ µ µ NZh∥X∥Zh∥X ,∀Zh∈Xh,(2.41) for a given UN∈XN. 56 2.3. A posteriori error estimator The construction of an a posteriori error estimator goes through the proof of the existence and uniqueness of the solution of problem (2.22) . Theorem 2.1 (Uniqueness).Let µ µ µ∈D , and assume that βN(µ µ µ)>0 . If problem (2.22) admits a solution Uh(µ µ µ)such that Uh(µ µ µ)∈BXUN(µ µ µ),βN(µ µ µ) ρT then, this solution is unique. For the proof of this theorem and the next one, we define the following operators: • The residual R(·;µ µ µ):Xh−→X′ hby ⟨R(Zh;µ µ µ),Vh⟩=A(Zh,Vh;µ µ µ)−F(Vh;µ µ µ),∀Zh,Vh∈Xh.(2.42) • The derivative of the operator A(·,·;µ µ µ),DA(Uh(µ µ µ);µ µ µ):Xh−→X′ h, by ⟨DA(Uh(µ µ µ);µ µ µ)Zh,Vh⟩=∂1A(Uh,Vh;µ µ µ)(Zh),∀Zh,Vh∈Xh.(2.43) • The operator H:Xh−→Xhby H(Zh;µ µ µ) = Zh−DA(UN(µ µ µ);µ µ µ)−1R(Zh;µ µ µ),∀Zh∈Xh.(2.44) Proof. The strategy of this proof is to show that H(·;µ µ µ) is a contraction in the sense of Definition A.3. If there exists a fixed point Uh , this point is a solution of problem (2.22) since H(Uh(µ µ µ);µ µ µ) = Uh(µ µ µ)−DA(UN(µ µ µ);µ µ µ)−1R(Uh(µ µ µ);µ µ µ) = Uh(µ µ µ) ⇐⇒R(Uh(µ µ µ);µ µ µ) = 0. Let us study if the operator H is a contraction, proving in this way that it has a unique fixed point. By definition, we have H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ)=(Z1 h−Z2 h)−DA(UN(µ µ µ);µ µ µ)−1(R(Z1 h;µ µ µ)−R(Z2 h;µ µ µ)).(2.45) Using the residual definition (2.42) and the mean value theorem, it holds R(Z1 h;µ µ µ)−R(Z2 h;µ µ µ) = DA(ξ;µ µ µ)(Z1 h−Z2 h),(2.46) 2. RB steady Smagorinsky model with variable geometry 57 where ξ=sZ1 h+(1−s)Z2 h, for some s∈(0,1). To prove this, we define T:[0,1]−→R, by T(s) = ⟨R(sZ1 h+(1−s)Z2 h;µ µ µ),Vh⟩,∀Vh∈Xh. The function Tis derivable for all s∈(0,1)and its derivative is T′(s) = ⟨DA(sZ1 h+(1−s)Z2 h;µ µ µ)(Z1 h−Z2 h),Vh⟩. If we apply the mean value theorem to Twe have (2.46). Now, multiplying (2.45) by DA(UN(µ µ µ);µ µ µ)and applying (2.46), we can write DA(UN(µ µ µ);µ µ µ)(H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ)) = [DA(UN(µ µ µ);µ µ µ)−DA(ξ;µ µ µ)](Z1 h−Z2 h). Then, thanks to the Lipschitz condition (2.34), we can write ⟨DA(UN(µ µ µ);µ µ µ)(H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ)),Vh⟩ ≤ρT∥UN(µ µ µ)−ξ∥X∥Z1 h−Z2 h∥X∥Vh∥X(2.47) Taking Zh=H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ)in (2.41), then, βN(µ µ µ)∥H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ)∥X∥Tµ µ µ N(H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ))∥X ≤⟨DA(UN(µ µ µ);µ µ µ)(H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ)),Tµ µ µ N(H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ))⟩, and applying (2.47) to Vh=Tµ µ µ N(H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ)) we have that, βN(µ µ µ)∥H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ)∥X∥Tµ µ µ N(H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ))∥X ≤ρT∥UN(µ µ µ)−ξ∥X∥Z1 h−Z2 h∥X∥Tµ µ µ N(H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ))∥X from which we obtain that ∥H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ)∥X≤ρT βN(µ µ µ)∥UN(µ µ µ)−ξ∥X∥Z1 h−Z2 h∥X. If Z1 h and Z2 h are in BX(UN(µ µ µ),α) then ∥UN(µ µ µ)−ξ∥X≤α since ξ=sZ1 h+(1−s)Z2 h for some s∈(0,1)and ∥H(Z1 h;µ µ µ)−H(Z2 h;µ µ µ)∥X≤ρT βN(µ µ µ)α∥Z1 h−Z2 h∥X Then, H(·;µ µ µ)is a contraction if α<βN(µ µ µ) ρT . In conclusion, if there exists a fixed point Uh in BXUN(µ µ µ),βN(µ µ µ) ρT , this point is unique. 64 2.4. Numerical results 0 10 20 30 40 50 60 70 M 10-4 10-3 10-2 10-1 100 max D|g( )-IM[g( )]| /|g( )| Figure 2.2: Convergence of the EIM algorithm. 10-12 3 10-10 2.9 10-8 4 10-6 2.8 3.5 10-4 2.7 3 10-2 2.6 2.5 2.5 2 2 4 6 8 10 12 10-4 Figure 2.3: Normalized error in L∞-norm of the turbulent eddy viscosity. 2.4.3 Setup for an a posteriori error bound estimate This section is intended to detail the computation of the different items that take part in the construction of the a posteriori error bound estimate defined in (2.50) . • For the computation of the Stability Factor (2.38) we substitute this quantity for 2. RB steady Smagorinsky model with variable geometry 65 βh(µ µ µ) since UN(µ µ µ) is intended to be a good approximation of Uh(µ µ µ) drawn from the solutions that we already computed for the EIM. To compute the inf-sup stability factor, we use the procedure exposed in Section 1.5 in [ 14 ] using the Radial Basis Function (RBF) algorithm to approximate the value of the inf-sup stability factor for the new parameters. First, we compute βh(µ µ µ) for the 81 solutions getting the lowest eigenvalue in problem (Find (α,Zh)∈R×Xh,Zh=0,such that F(µ)TX−1F(µ)Zh=αXZh,∀Zh∈Xh, where the matrices X and F(µ) are the matrices associated to the inner product related to the X-norm (2.13) and the tangent operator ∂1A, this is, VT hXZh= (Vh,Zh)X,∀Vh,Zh∈Xh, VT hF(µ)Zh=∂1A(Uh(µ),Vh;µ)(Zh),∀Vh,Zh∈Xh,∀µ∈D, resulting that βh(µ µ µ) = (αmin)1/2. Secondly, we apply the RBF algorithm to obtain an approximation of βh(µ µ µ) for all µ µ µ∈D . In this case, we stop the algorithm when the estimator is below εβ=10−4 . In the RBF, we tested 729 solutions and at the end we have selected 82 parameters. • The residual (2.49) is computed in each iteration solving a variational problem applying the Riesz representation Theorem A.5. • To compute the Lipschitz constant ρT defined in (2.34) we need to compute η∗ , λ∗ , ϕ∗ , ˆ ϕ and h∗ max which are easy to obtain since these parameters only depend on D and C4;T defined in (2.30) . This last constant is approximated by a fixed-point algorithm described in [15, 34]. For the first iterations, τN(µ µ µ) is not lower than one since the residual is still large. While τN(µ µ µ)>1 , we use as a posteriori error bound estimator the proper τN(µ µ µ) . 2.4.4 Offline phase Finally, we compute the RB space applying the Greedy algorithm for the selection of dim(Dtrain) = 625 . We stop the algorithm when ∆N(µ µ µ)<εRB with εRB =5·10−4 . In each iteration, we use the Gram-Schmidt algorithm to orthonormalize the reduced basis functions to improve the condition number of the reduced problem. 66 2.4. Numerical results 5 10 15 20 25 N 10-4 10-3 10-2 10-1 100 101 102 103 104 max D N( ) max D N( ) Figure 2.4: Greedy Algorithm convergence. Finally, in Figure 2.4 we show the Greedy Algorithm convergence. We observe that τN(µ µ µ) is lower than one when we achieve 25 basis functions, then we are able to compute ∆N(µ µ µ) finally stopping the algorithm with 28 basis function. It is interesting to take a look at the parameter selection. In Figure 2.5 we show the a posteriori error bound estimator for six different iterations of the Greedy Algorithm. As we have explained before, in figures 2.5a-2.5d we show τN(µ µ µ) while in figures 2.5e and 2.5f we show ∆N(µ µ µ) . We observe the progressive decrease of errors. The Greedy algorithm starts with µ µ µ0= (ωmin,σmin) as the first parameter, and in the following iterations it selects µ µ µ1= (ωmax,σmax) , µ µ µ2= (ωmax,σmin) and µ µ µ3= (ωmin,σmax) finishing with µ µ µ4= (2.771,3) that is close to (ωmean,σmean) . This situation is expected since the algorithm is catching the main information in these iterations. This is described in Figure 2.5a. In Figure 2.5b we still see a homogeneous distribution of the parameter selection and that the total error is decreasing along figures 2.5c and 2.5d. Finally, in Figure 2.5e we use the a posteriori error bound estimator and we see that in three more iterations, the maximum is lower than εRB and we finish the algorithm (Figure 2.5f). 2. RB steady Smagorinsky model with variable geometry 67 10-10 3 10-5 2.9 4 2.8 100 3.5 2.7 3 105 2.6 2.5 2.5 2 50 100 150 200 250 (a) Iteration 5. 10-10 3 10-5 2.9 4 2.8 100 3.5 2.7 3 105 2.6 2.5 2.5 2 10 20 30 40 50 60 (b) Iteration 10. 10-10 3 10-5 2.9 4 2.8 100 3.5 2.7 3 105 2.6 2.5 2.5 2 1 2 3 4 5 6 7 (c) Iteration 15. 10-10 3 10-5 2.9 4 2.8 100 3.5 2.7 3 105 2.6 2.5 2.5 2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 (d) Iteration 20. 10-10 3 10-8 2.9 4 10-6 2.8 3.5 10-4 2.7 3 10-2 2.6 2.5 2.5 2 1 2 3 4 5 6 7 8 9 10-4 (e) Iteration 25. 10-10 3 10-8 2.9 4 10-6 2.8 3.5 10-4 2.7 3 10-2 2.6 2.5 2.5 2 0.5 1 1.5 2 2.5 3 3.5 4 4.5 10-4 (f) Iteration 28. Figure 2.5: ∆N(µ µ µ)through the Greedy Algorithm. 68 2.4. Numerical results Data Case 1 Case 2 Case 3 ∆N1.21·10−41.07·10−41.89·10−4 ∥Uh−UN∥X5.93·10−63.73·10−67.28·10−6 Efficiency 20.4 28.69 25.96 speedup 142 167 148 Table 2.1: Errors and speedups. We have tried to keep going with the algorithm looking for decreasing the error to εRB =10−4 , however, in a few more iterations, the algorithm has reached a point when it selects parameters that already have been selected. This makes the reduced system be singular, thus stopping the algorithm. This can be due to a low precision of the solution of the linear problems that appear in the full order model. 2.4.5 Online phase In this section, we validate an a posteriori error bound estimator ∆N(µ µ µ) and we compare the computational time between the Reduced Basis and the Finite Element solutions. To do this, we choose three random pairs of parameters: • Case 1: ω=2.891m, σ=2.734m, • Case 2: ω=2.649m, σ=2.65m, • Case 3: ω=2.469m, σ=2.923m, and we compute the Uh(µ µ µ) and UN(µ µ µ) for each one studying the error between them. This computation has been made on a Mac Book Pro 2017 with a 2.3 GHz Intel Core i5 processor. The computational time for the Finite Element (FE) solutions Uh(µ µ µ) is between 3 and 3.60min while for the Reduced Basis solutions UN(µ µ µ) is between 1.3 and 1.48s , which is a considerable improvement. In the Table 2.1 we show the value of the estimator ∆N(µ µ µ) , the error between both solutions in the X -norm, the estimator efficiency and the computing time relation between the FE and the RB solutions. The speed-ups factors are close to 150 . Finally, in Figures 2.6-2.7 we show a comparison between the estimate ∆N(µ µ µ) and the real error ∥Uh−UN∥X . As we can see, both have the same shape and the real error is lower than the estimate, which validates the estimate and the computation made in this section. Additionally, in figures 2.8 and 2.9 we show the velocity field and the pressure, respectively, for the three cases. The domain represented in these figures is Ω(µ µ µ) for each case. 2. RB steady Smagorinsky model with variable geometry 69 10-10 3 10-8 4 2.8 10-6 3.5 Estime for velocity and pressure 10-4 3 2.6 2.5 2.4 2 0.5 1 1.5 2 2.5 3 3.5 4 4.5 10-4 Figure 2.6: Final estimate ∆N(µ µ µ). 10-10 3 10-8 4 2.8 10-6 3.5 Error between FE and RB for vel and pres in T-norm 10-4 3 2.6 2.5 2.4 2 2 4 6 8 10 12 14 16 18 10-6 Figure 2.7: Final error ∥Uh−UN∥X. 70 2.4. Numerical results Figure 2.8: Velocity field for each case. Figure 2.9: Pressure for each case. 3 Reduced Basis steady Smagorinsky model for forced convection flow In this chapter, we afford the thermal comfort optimization of a cloister by means of Reduced Basis modeling. We first build the RB model and then apply it to the thermal comfort optimization of the cloister. We consider forced convection, in which the momentum conservation equations are decoupled from the temperature equation. We introduce the new problem in Section 3.1, bringing in the temperature to the Smagorinsky model with the Boussinesq approximation, defining the related RB problem, and explaining the connection with the results presented in Chapter 2. Next, we build an a posteriori error bound estimator for the temperature in Section 3.2, following an classical technique using the Lax-Milgram Lemma A.1. Finally, we show some numerical results in Section 3.3 for a realistic case validating the estimate. In Section 3.4, we also solve an optimization problem which goal is to choose the pair of parameters that give us the best cloister configuration minimizing the minimizing the difference between the temperature and an ideal comfort temperature near the ground. 3.1 Problem statement In this section, we introduce an aero-thermal model for the cloister. We modify the Smagorinsky model (2.1) including the Boussinesq approximation for temperature. 71 72 3.1. Problem statement We consider forced convection in which buoyancy forces are not considered in the momentum conservation equation. This is an acceptable approximation for moderate temperature gradients between the ground and the outer air. We refer to Chapter 14 in [ 49 ] for more information. In principle, a steady thermal flow governed by the Smagorinsky turbulence model would obey to the following equations:                                      −ν∆w+(w·∇)w−∇·(νt(w;µ µ µ)∇w)+∇p=−gβ(θ−θ0),in Ω(µ µ µ), ∇·w=0,in Ω(µ µ µ), (w·∇)θ−∇·((κ+κt(w))∇θ) = Q,in Ω(µ µ µ), w= (g(·;µ µ µ),0),on ΓS(µ µ µ), w=0,on ΓW(µ µ µ)∪ΓF(µ µ µ), θ=θS,on ΓS(µ µ µ), θ=θF,on ΓF(µ µ µ), −n·((κ+κt(w;µ µ µ))∇θ) = α ρcp (θ−θ0),on ΓW(µ µ µ). (3.1) An additional term appears in the momentum conservation equation (3.1)1 , −gβ(θ−θ0) , that models the buoyancy force, where β is the coefficient of expansion of the fluid, g is the gravitational acceleration, and θ0 represent a reference temperature. In the energy conservation equation (3.1)3 , the first term represents the transport of heat due to the motion of the fluid and it is known as the advection term, while the second one is the diffusion term where we introduce the eddy diffusivity κt(·;µ µ µ) , which usually in turbulence modeling is assumed to be proportional to νt(·;µ µ µ) . At the right-hand side we find the heat generation term Q . The coefficient κ is known as the thermal diffusivity and is defined as κ=k/ρcp , where k is thermal conductivity, ρ is mass density and cp is the specific heat capacity (see Table 1). Finally, we impose Dirichlet boundary conditions at the top of the domain ( ΓS(µ µ µ) ) and on the floor ( ΓF(µ µ µ) ). Moreover, we impose the heat flow on the walls ( ΓW(µ µ µ) ), where α is the heat transfer or film coefficient. This is stated in the last equations (3.1)6 - (3.1)8 . For simplicity, we do not consider any distributed heating for this problem, then Q=0 . For a perfect gas, β=1/θ , then the buoyancy term is written as −g(θ−θ0)/θ . If the absolute variation of the temperature from the ground to the outer air is sufficiently small, then the buoyancy effects can be neglected, and we obtain a variation of the temperature depending on the velocity, not the other way round. This is the situation that we consider, called “forced convection”. 3. RB steady Smagorinsky model for forced convection flow 73 To sum up, we assume forced convection, which gives us the following model (1)           −ν∆w+(w·∇)w−∇·(νt(w;µ µ µ)∇w)+∇p=0,in Ω(µ µ µ), ∇·w=0,in Ω(µ µ µ), w= (g(·;µ µ µ),0),on ΓS(µ µ µ), w=0,on ΓW(µ µ µ)∪ΓF(µ µ µ), (3.2) (2)                (w·∇)θ−∇·((κ+κt(w;µ µ µ))∇θ) = 0,in Ω(µ µ µ), θ=θS,on ΓS(µ µ µ), θ=θF,on ΓF(µ µ µ), −n·((κ+κt(w;µ µ µ))∇θ) = α ρcp (θ−θ0),on ΓW(µ µ µ). (3.3) where we first compute the velocity w and the pressure p in (3.2) and in a second step, we compute the temperature θ in (3.3) . As the velocity does not depend on the temperature, we just need to build a RB problem for the temperature. We start by obtaining the dimensionless problem following the same procedure as in Section 2.1. We consider θ0=θC where θC is an ideal comfort temperature and we set the dimensionless temperature θ∗=θ/θC . We obtain the dimensionless equations                          (w∗·∇∗)θ∗−∇∗·1 Pe +κ∗ t(w∗;µ µ µ∗)∇∗θ∗=0,in Ω(µ µ µ∗), θ∗=θS θC ,on ΓS(µ µ µ∗), θ∗=θF θC ,on ΓF(µ µ µ∗), −n·1 Pe +κ∗ t(w∗;µ µ µ∗)∇∗θ∗=Nu Pe (θ∗−1),on ΓW(µ µ µ∗), (3.4) where w=U0w∗ , κ∗ t(w∗;µ µ µ∗) = 1 Prν∗ t(w∗;µ µ µ∗), where we recall that ν∗ t(w∗;µ µ µ∗) was defined in (2.2) , Pr , Pe and Nu are the Prandlt, Péclet and Nusselt numbers defined in Appendix (A.21) , (A.23) and (A.24) respectively, such that Pr =cpµ k,Pe =RePr,Nu =αL k, where µ is the dynamical viscosity. Again, from now on, we avoid the star notation for simplicity. 80 3.2. A posteriori error estimator and 4 ∑ l=1 2 ∑ k=1 γl kkaθ t,kl(wN,θuN h,θuN h;µ µ µ)≥ˆ γ∑ K∈ThZK νt(wN;µ µ µ) Pr |∇θuN h|dK ≥ˆ γ Pr min x∈Ωνt(wN;µ µ µ)∥∇θuN h∥2 L2(Ω) We bound the boundary term by 4 ∑ l=1|JΦl|dθ l(θuN h,θuN h)≥min l=1,3,4|JΦl|ZΓW|θuN h|2dΓ=min l=1,3,4|JΦl|∥θuN h∥2 L2(ΓW)≥0. Finally, we upper bound the last term and we prove that it is sufficiently small not to affect the coercivity. We have 4 ∑ l=1 2 ∑ k=1 ηl kkcθ kl(wN,θuN h,θuN h) = ZΩ(µ µ µ)(wN·∇)(θuN h)2dΩ =−1 2ZΩ(µ µ µ) ∇·wN(θuN h)2dΩ(3.14) Since wN=uN+uD and uN is the solution of problem (2.26) , it verifies that for all qN∈QN,ZΩ(µ µ µ) ∇·wNqN=0. Therefore, wN is weakly divergence free and we can not ensure that (3.14) is null, since we can not guarantee that (θuN h)2 belongs to QN . However, w solution of (2.1) is divergence free, then ZΩ(µ µ µ) ∇·wN(θuN N)2dΩ =ZΩ(µ µ µ) ∇·(wN−w)(θuN N)2dΩ ≤∥∇(wN−w)∥L2(Ω(µ µ µ))∥(θuN N)2∥L2(Ω(µ µ µ)) =∥wN−w∥H1 0(Ω(µ µ µ))∥θuN N∥2 L4(Ω(µ µ µ)) Hence, ∥θuN N∥2 L4(Ω(µ µ µ)) =ZΩ(µ µ µ)(θuN N)41/2 = 4 ∑ l=1ZΩl|JΦl|(θuN N)4!1/2 ≤max l=1,2,3,4|JΦl|1/2ZΩ(θuN N)41/2 =max l=1,2,3,4|JΦl|1/2∥θuN N∥2 L4(Ω) ≤C2 4;1,2max l=1,2,3,4|JΦl|1/2∥θuN N∥2 H1 0(Ω) 3. RB steady Smagorinsky model for forced convection flow 81 where the constant C4;1,2is defined in (3.11). Then, 4 ∑ l=1 2 ∑ k=1 ηl kkcθ kl(wN,θuN h,θuN h)>−εN(µ µ µ) 2∥θuN N∥2 H1 0(Ω) where εN(µ µ µ) = C2 4;1,2max l=1,2,3,4|JΦl|1/2∥wN−w∥H1 0(Ω(µ µ µ)). Finally, the coercivity constant is β(µ µ µ) = ˆ γ1 Pe +min x∈Ω νt(wN;µ µ µ) Pr −εN(µ µ µ) 2, which is positive for small enough ∥wN−w∥H1 0(Ω(µ µ µ)). Remark 3.2. If N is sufficiently large, then we can obtain that ∥wN−w∥H1 0(Ω(µ µ µ)) to be sufficiently small to guarantee that β(µ µ µ) be positive and in this way, to ensure the coercivity of the form Aθ. Once we obtain the continuity and coercivity factors, we are able to infer an a posteriori error estimator. Lemma 3.1. Let µ µ µ∈D . We denote euN h(µ µ µ) = θuN h(µ µ µ)−θuN N(µ µ µ)∈Θh the error between the high fidelity solution of problem (3.9) and the reduced solution of (3.10). Then, ∥R(wN;θuN N(µ µ µ))∥Θ′ h γh(µ µ µ)≤∥euN h(µ µ µ)∥H1 0(Ω)≤∥R(wN;θuN N(µ µ µ))∥Θ′ h β(µ µ µ),(3.15) where the residual operator R(wN;θuN N(µ µ µ)) ∈Θ′ his defined by ⟨R(wN,θuN N(µ µ µ));θv h⟩=Fθ(wN,θv h;µ µ µ)−Aθ(wN,θuN N(µ µ µ),θv h;µ µ µ),∀θv h∈Θh. Proof. We represent R(wN;·)in Θhby means of the L2-inner product by r(wN)∈Θh: (r(wN),θv h)Ω=⟨R(wN;θuN N(µ µ µ)),θv h⟩,∀θv h∈Θh.(3.16) Now taking into account that the problem defined by euN hcan be expressed by Aθ(wN,euN h,θv h;µ µ µ) = Aθ(wN,θuN h,θv h;µ µ µ)−Aθ(wN,θuN N,θv h;µ µ µ) =Fθ(wN,θv h;µ µ µ)−Aθ(wN,θuN N,θv h;µ µ µ) =⟨R(wN;θuN N(µ µ µ)),θv h⟩. (3.17) 82 3.3. Numerical results and thanks to the continuity of the operator Aθ(wN,·,·;µ µ µ)in Theorem 3.1, it holds |(r(wN),θv h)Ω|=|Aθ(wN,euN h,θv h;µ µ µ)|≤γh(µ µ µ)∥euN h∥H1 0(Ω)∥θv h∥H1 0(Ω). Due to the definition of dual norm, we obtain that ∥R(wN;θuN N(µ µ µ))∥Θ′ h≤γh(µ µ µ)∥euN h∥H1 0(Ω). On the other hand, we apply the Lax-Milgram Lemma A.1 to the problem (3.17), ∥euN h∥H1 0(Ω)≤∥R(wN;θuN N(µ µ µ)∥Θ′ h β(µ µ µ). Remark 3.3. In view of the preceding lemma, we state an a posteriori error estimator as follows ∆θ,N(µ µ µ) = ∥R(wN;θuN N(µ µ µ)∥Θ′ h β(µ µ µ).(3.18) Then the estimation (3.15) can be rewritten as ∆θ,N(µ µ µ)β(µ µ µ) γh(µ µ µ)≤∥euN h(µ µ µ)∥H1 0(Ω)≤∆θ,N(µ µ µ). The efficiency is the quotient between the upper and lower bounds of the error, and in this case is γh(µ µ µ)/β(µ µ µ)for all µ µ µ∈D. 3.3 Numerical results In this section, we apply the RB method to problem (3.9) computing the RB spaces using the estimate (3.18) . We consider the velocity field uN∈YN solution of (2.26) , for this reason we keep the same domain configuration as in Section 2.4.1. Again, we recall that wN=uN+uD . 3.3.1 Problem statement We set the floor temperature equal to the comfort one that we assume to be 295.15K ( 22°C ). We assume that the temperature of the outer air is 303.15K ( 30°C ). Since the characteristic temperature was chosen as θC=θF , the dimensionless temperature on the floor should be 1 , and 1.4 on the top. We set the Prandlt number associated to the air equal to 0.71 and the Nusselt number associated to the transfer of temperature at the walls equal to 2 (cf. [ 4 ]). The Péclet number appearing in the temperature equation is defined by Pe =RePr . Considering that Re =3100 , then Pe =2201 . Finally, we set the parameter set D as in Section 2.4.1, that is µ µ µ= (ω,σ)∈D= [2m,4m]×[2.5m,3m] . 3. RB steady Smagorinsky model for forced convection flow 83 3.3.2 Offline phase We again use the EIM applied to approximate the eddy diffusion, taking into account the approximation of the eddy viscosity, already computed in Chapter 2. In this case, the computation of an a posteriori estimator is simpler than the one developed in Chapter 2. Besides, we need to compute an approximation for the coercivity constant for each µ µ µ∈D and the residual associated to the problem (3.9) . To compute the coercivity constant, we apply the same procedure as for the inf-sup constant in Chapter 2. We compute βh(µ µ µ) for 81 pairs of parameters with the procedure described in Section 3.7 in [ 38 ], solving an eigenvalue problem, and then, we apply the RBF algorithm explained in Section 1.5 in [ 14 ]. We stop the algorithm when the estimator is below εβ=10−4 testing 729 parameter pairs and selecting 82 parameters. For the residual, we use a representation r(uN) in Θh as in (3.16) , solving a FE problem. Finally, we compute the offline phase stopping the algorithm when ∆θ,N(µ µ µ)<εθ,RB with εθ,RB =5·10−3 . Again, we orthonormalize the RB functions using the Gram-Schmidt algorithm with H1 0(Ω) -norm to the temperature space Θh . Along this phase, we have tested a total of 625 parameters. In Figure 3.1 we can see the Greedy Algorithm convergence. It stops when 17 basis functions are selected. 3.3.3 Online Phase In this section, we compare the computational time for the full discrete model and the RB model, and we validate the a posteriori error bound estimator ∆θ,N(µ µ µ) . We select the three random cases, already chosen in Section 2.4.5: • Case 1: ω=2.891m, σ=2.734m, • Case 2: ω=2.649m, σ=2.65m, • Case 3: ω=2.469m, σ=2.923m. First, we compute the high fidelity problem (computing θuh h(µ µ µ) from uh(µ µ µ) solution of (2.22) ) and the full reduced problem (computing θuN N(µ µ µ) from uN(µ µ µ) solution of (2.26) ). We measure the error between both solutions and the time that takes to compute them and we compare them to obtain the speedup. In the Table 3.1 we can observe the committed error in H1 0(Ω) -norm and the speedup. As we can see, the errors are admissible and we reduce the computational time from 84 3.3. Numerical results 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 N 10-3 10-2 10-1 100 101 102 max D ,N( ) Figure 3.1: Greedy Algorithm convergence. Data Case 1 Case 2 Case 3 ∥θuh h−θuN N∥H1 0(Ω)5.10403·10−46.24142·10−41.09348·10−3 speedup 134 146 156 Table 3.1: Errors and speedups. 3.5min to 1.33s . This computation has been done on only one processor of a cluster with CPUs AMD EPYC 7542 2.9 GHz. In order to validate the a posteriori error bound estimator ∆θ,N(µ µ µ) , we consider the finite element solution θuN h(µ µ µ) computed from the velocity uN(µ µ µ) and we compare the error between θuN h(µ µ µ) and θuN N(µ µ µ) in the H1 0(Ω) -norm and the estimator ∆θ,N(µ µ µ) for 625 pair of parameters. In Figure 3.2 we show two graphics in 3D with the same limits in the axes. In Figure 3.2a the vertical scale goes from 5.03 ·10−9 to 5·10−3 and in Figure 3.2b it goes from 1.61·10−9 to 1.6·10−3 . Both graphics have the same shape and the estimate is above the error between the finite element and the RB solution. With this, we can validate our estimate for the temperature. Additionally, we compute the estimate and the error for the three previous cases of parameters. The full high fidelity model computes the temperature associated to uh(µ µ µ) , 3. RB steady Smagorinsky model for forced convection flow 85 10-8 3 10-6 4 2.8 10-4 Estime for temperature 3.5 3 2.6 2.5 2.4 2 0.5 1 1.5 2 2.5 3 3.5 4 4.5 10-3 (a) Final ∆θ,N(µ µ µ). 10-8 3 10-6 4 2.8 10-4 Exact error between FE and RB for temp 3.5 3 2.6 2.5 2.4 2 2 4 6 8 10 12 14 10-4 (b) ∥θuN h−θuN N∥H1 0(Ω). Figure 3.2: Validation of the estimate. that is, θuh h(µ µ µ) , however, we are using the temperature associated to uN(µ µ µ) , that is, θuN h(µ µ µ) as the high fidelity solution. We show the results in Table 3.2 and we also see that ∥θuh h−θuN N∥H1 0(Ω) is close to ∥θuN h−θuN N∥H1 0(Ω) at least for these three cases. Data Case 1 Case 2 Case 3 ∥θuN h−θuN N∥H1 0(Ω)5.10131·10−46.30339·10−41.11298·10−3 ∆θ,N1.76917·10−32.78895·10−33.81012·10−3 Efficiency 3.47 4.42 3.42 Table 3.2: Validation for the three cases of the estimate. Finally, in Figure 3.3 we show the temperature in Kelvin ( K ) for the three cases. It is interesting to compare them with the velocity field shown in Figure 3.4. We observe that the highest temperate values are in the biggest eddy, while the temperature keeps colder in the corridors as expected from buoyancy effects, even if these have not been taken into account to compute the velocity. 3.4 Application: geometrical optimization of a cloister In this section, we optimize the cloister geometry with the goal of maximizing the comfort measured as the difference between an ideal comfort temperature and the computed temperature, in a region where the visitors of the cloister may find themselves. To do 86 3.4. Application: geometrical optimization of a cloister Figure 3.3: Temperature for the three test cases. Figure 3.4: Velocity for the three test cases. 3. RB steady Smagorinsky model for forced convection flow 87 this, we define the “thermal comfort” functional J:D−→R+ by J(µ µ µ) = 1 A(ΩI(µ µ µ)) ZΩI(µ µ µ)|θ(µ µ µ)−θC|2dΩ(3.19) where A(Ω) represents the area of Ω , ΩI(µ µ µ) = Sl∈IΩl⊂Ω(µ µ µ) , θ(µ µ µ) is the solution of (3.3) and I represents the set of indices of the domains where the visitors of the cloister may find them selves. The purpose is to minimize this functional for µ µ µ∈D and choosing ΩI depending on our interest: • we choose I={1,2,3}if we are interested on the ground floor, • we choose I={1,3}if we are interested on the corridors. We have to solve the constrained minimization problem µ µ µopt =arg min µ µ µ∈Dad JN(µ µ µ) where Dad is the set of admissible parameters from the architectural point of view. Instead of applying optimization methods to solve the optimization problem, we compute the functional for a training set Dtrain ⊂Dad and we directly look for the minimum. This is possible thanks to the RB solution of the temperature equation, whose computational cost is very low. With this method, we obtain additional information about the functional based on the parameters and in particular, the minimum. Since we assume that θuN h(µ µ µ) is a good approximation of θ(µ µ µ) solution of (3.3) and we build θuN N(µ µ µ) from θuN h(µ µ µ) , we can say that we are not loosing the important information on the temperature behavior using the RB method for (3.19) with θ=θuN N , denoted by JN(µ µ µ) . Again, we select a total of 625 pairs of parameters and we compute JN(µ µ µ) , this takes less than 16 minutes of computing time. We remember that θuN N(µ µ µ) solution of (3.10) is a dimensionless variable defined on Ωr with θuN N=0 on ΓS∪ΓF , then we need to set the Dirichlet boundary condition θD to θuN N(µ µ µ) . In practice, we need to define a dimensionless functional J∗ N(µ µ µ) and for simplicity, we define it on ΩI(µ µ µr) . On one side, we know that the physical dimensions of the functional in (3.19) is θ2 0 , then we define J∗ N=JN/θ2 C . On the other side, we use the mapping defined in (2.5) to define the functional on ΩI(µ µ µr) . Finally, we obtain that J∗ N(µ µ µ) = 1 A∗(I;µ µ µ)∑ i∈IZΩi|JΦi||θuN N(µ µ µ)+θD−1|2dΩ,(3.20) 88 3.4. Application: geometrical optimization of a cloister 2050 4 2100 2150 2200 3.5 3 2250 2.9 2300 3 2350 2.8 2400 2.7 2.5 2.6 22.5 2100 2150 2200 2250 2300 2350 K (a) JN(µ µ µ)for I={1,2,3}. 150 4 200 250 300 3.5 3 350 400 2.9 450 3 500 2.8 550 2.7 2.5 2.6 22.5 200 250 300 350 400 450 500 K (b) JN(µ µ µ)for I={1,3}. Figure 3.5: Computation of the functional for θuN N(µ µ µ). where A∗(I;µ µ µ) = ∑ i∈IZΩi|JΦi|dΩ=(σWif I={1,2,3} 2σω if I={1,3}.(3.21) In Figure 3.5 we show the functional JN(µ µ µ) for both cases. In these graphics, we can see that the functional minimum it is achieved for a corridor width of 4m and a height of 2.5m , which correspond to the maximum width and the minimum height of the corridor. This makes sense since, with this configuration, the air barely flows through the corridors, it stays in the interior of the courtyard and does not carry the hot air that comes from the top of the corridors. In the hypothetical case that the corridors were higher, the horizontal flow that comes into the corridors would be increased, and so would be the temperature. Similarly, if the corridors were narrower, their area would decrease and so would the amount of air in them, on the contrary, it would not affect the flow that comes into the corridors. A Mac Book Pro 2017 with a 2.3 GHz Intel Core i5 processor takes 7.79min for case I={1,2,3} and 8.04min for case I={1,3} using two processors. Open problems To consider more realistic situations, we set some open problems related to the modeling carried on this part: □ Extension to free convection, where buoyancy effects are included in the momentum conservation equation. An estimator involving velocity, pressure, and temperature must be developed. □ Cloister modeling in 3D. It may be necessary to change the mapping Φ developed in (2.5) and the corresponding a posteriori error estimate. □ Inclusion of radiation effects Q to the problem using the previous modeling in 3D. This is necessary since radiation depends on the position of the sun and this movement is in 3D. □Inclusion of unsteady and seasonal effects. □Inclusion of tempering elements in the cloister, such as fountains, vegetation, etc. 96 4.1. Problem statement and the nonlinear Smagorinsky term, at(·;·,·) is defined by at(u;z,v) = ZΩ νt(u)∇z:∇vdΩ.(4.8) To simplify the notation, we shall consider the following discrete functions: •˜ vh:(0,Tf)−→Yhis the piecewise constant function defined by ˜ vh(t):=vk hif t∈(tk−1,tk). • ˜qh:(0,Tf)−→Qhis the piecewise constant in time function defined by ˜qh(t):=qk hif t∈(tk−1,tk). •∂⋆ t˜ vh:(0,Tf)−→Yhis the discrete derivative of ˜ vhdefined by ∂⋆ t˜ vh(t):=vk h−vk−1 h ∆tif t∈[tk−1,tk]. •˜ f:(0,Tf)−→L2(Ω) is the piecewise constant in time function defined by ˜ f(t) = fk if t∈(tk−1,tk). We can define ˜ Yh as the space formed by the piecewise constant functions that belongs to Yh and ˜ Qh formed by the piecewise constant functions that belongs to Qh . Furthermore, we can define an equivalent problem to (4.4) as follows: Find (˜ uh,˜ph)∈˜ Yhט Qh solution of          ˜m(∂∗ t˜ uh,˜ vh)+ 1 µ˜a(˜ uh,˜ vh)+ ˜at(˜ uh;˜ uh,˜ vh) +˜c(˜ uh,˜ uh,˜ vh)+ ˜ b(˜ vh,˜ph) = ⟨˜ f,˜ vh⟩QT,∀˜ vh∈˜ Yh, ˜ b(˜ uh,˜qh) = 0,∀˜qh∈˜ Qh, (4.9) where the forms ˜m(·,·) , ˜a(·,·) , ˜ b(·,·) , ˜c(·,·,·) and ˜at(·;·,·) are the same as the forms defined in (4.5) - (4.8) by changing the integration domain by the time-space domain QT=If×Ω , instead of Ω . Note that for a given Banach Space X and any n∈N , the continuous norms applied to the piecewise functions ˜ v:(0,Tf)−→X such that ˜ v(t):=vk , for all t∈(tk−1,tk) and k=1,...,L are related to its discrete norms, this is, ∥˜ v∥Ln(X)=ZTf 0∥˜ v(s)∥n Xds1/n = L ∑ k=1Ztk tk−1∥˜ v(s)∥n Xds!1/n = L ∑ k=1 ∆t∥vk∥n X!1/n , and ∥˜ v∥L∞(X)=max t∈(0,Tf)∥˜ v(t)∥X=max k=1,...,L∥vk∥X. 4. Numerical approximation of transient Smagorinsky model 97 Remark 4.2. We shall denote by C a constant that may change from one to other occurrence, but which always is independent of h. Remark 4.3. Note that, since fk=1 ∆tZtk tk−1 f(s)ds for k=1,...,L , then ∥˜ f∥Ln(X)≤∥f∥Ln(X) for any Banach space X . This result is well know, however, we follow with the proof for more completeness. This is, for any n ≥1 ∥˜ f∥n Ln(X)=ZTf 0∥˜ f(s)∥n Xds = L ∑ k=1Ztk tk−1∥˜ f(s)∥n Xds = L ∑ k=1 ∆t∥fk∥n X = L ∑ k=1 ∆t1−n   Ztk tk−1 f(s)ds    n X≤ L ∑ k=1 ∆t1−nZtk tk−1∥f(s)∥Xdsn . Using the Hölder’s Inequality A.1 on the interval (tk−1,tk), ∥˜ f∥n Ln(X)≤ L ∑ k=1 ∆t1−nZtk tk−1 1dsn−1Ztk tk−1∥f(s)∥n Xds= L ∑ k=1Ztk tk−1∥f(s)∥n Xds =∥f∥n Ln(X). 4.2 A priori error estimates for the velocity In this section, we will derive some a priori error estimates for the velocity in terms of norms in time and space. Lemma 4.1. Assume that f∈L1(0,Tf;L2(Ω)) and (˜ uh,˜ph)∈˜ Yhט Qh is a solution of (4.9). Then this solution satisfies ∥˜ uh∥L∞(L2)≤C∥f∥L1(L2).(4.10) Proof. For any k=1,...,L , (uk h,pk h)∈Yh×Qh is a solution of (4.4) . We take vk h=uk h and we obtain from the first equation in (4.4), 1 ∆tm(uk h,uk h)−1 ∆tm(uk−1 h,uk h)+ 1 µa(uk h,uk h) +at(uk h;uk h,uk h)+b(uk h,pk h)+c(uk h,uk h,uk h) = ⟨fk,uk h⟩Ω We note that a(uk h,uk h) = ∥∇uk h∥2 L2(Ω),c(uk h,uk h,uk h) = 0, and that b(uk h,pk h) = 0 since uk hverifies the second equation in (4.4). 98 4.2. A priori error estimates for the velocity Furthermore, using the property (A.3) and multiplying by 2∆twe obtain: ∥uk h∥2 L2(Ω)−∥uk−1 h∥2 L2(Ω)+∥uk h−uk−1 h∥2 L2(Ω)+21 µ∆t∥∇uk h∥2 L2(Ω) +2∆tat(uk h;uk h,uk h) = 2∆t⟨fk,uk h⟩Ω(4.11) Then, as at(uk h;uk h,uk h)≥0, we have ∥uk h∥2 L2(Ω)−∥uk−1 h∥2 L2(Ω)≤2∆t⟨fk,uk h⟩Ω Summing for k=1,...,mfor some m≥1, we obtain: ∥um h∥2 L2(Ω)≤2 m ∑ k=1 ∆t⟨fk,uk h⟩Ω≤2 m ∑ k=1 ∆t∥fk∥L2(Ω)∥uk h∥L2(Ω)≤2 m ∑ k=1 ∆t∥fk∥L2(Ω)∥˜ uh∥L∞(L2) Taking the maximum in m=1,...,L: ∥˜ uh∥L∞(L2)≤2 L ∑ k=1 ∆t∥fk∥L2(Ω)=2∥˜ f∥L1(L2)≤2∥f∥L1(L2) and therefore we obtain the result for C=2. To take into account the eddy diffusion effects in our estimates, we consider the following norm. Let X⊆H1 0(Ω) be a Banach Space, we define the turbulence inner product (·,·)T in this space as (v,z)T=ZΩ1 ¯ µ+ν⋆ t∇v:∇zdΩ,∀v,z∈X(4.12) where ν⋆ t=νt(u(¯ µ)) , with ¯ µ=arg min µ∈D(∑ K∈Th (CShK)2min x∈K|∇u(x,µ)|K|XK), where u(µ) is the minimum on time of the solution u(·,µ) of (4.3) for µ∈D . Remark 4.4. The norm associated to the inner product defined in (4.12) is equivalent to the H1 0(Ω)-norm. We shall denote by ∥· ∥Lp T(H1 0) the norm in Lp(0,Tf;H1 0(Ω)) defined by ∥v∥Lp T(H1 0)=ZTf 0∥v(t)∥p Tdt1/p ,for 1 ≤p<+∞,(4.13) ∥v∥L∞ T(H1 0)=max t∈(0,Tf)∥v(t)∥T,(4.14) where the T -norm is induced by the inner product (4.12) . 4. Numerical approximation of transient Smagorinsky model 99 Lemma 4.2. Under the same assumptions of Lemma 4.1, it holds ∥˜ uh∥2 L∞(L2)+1 µ∥∇˜ uh∥2 L2(QT)+(CShmin)2∥∇˜ uh∥3 L3(QT)≤C∥f∥2 L1(L2).(4.15) where hmin is the minimum size of all elements in the mesh Th. Proof. From (4.11), we sum for k=1,...,m: ∥um h∥2 L2(Ω)+ m ∑ k=1 2∆tµ∥∇uk h∥2 L2(Ω)+at(uk h;uk h,uk h)≤ m ∑ k=1 2∆t⟨fk,uk h⟩Ω We observe that at(uk h;uk h,uh) = ∑ K∈Th (CShK)2ZK|∇uk h|3dK ≥(CShmin)2∥∇uk h∥3 L3(Ω). In particular, ∥um h∥2 L2(Ω)+ m ∑ k=1 ∆t µ∥∇uk h∥2 L2(Ω)+ m ∑ k=1 ∆t(CShmin)2∥∇uk h∥3 L3(Ω)≤ m ∑ k=1 2∆t⟨fk,uk h⟩Ω. Moreover, thanks to (4.10) m ∑ k=1 ∆t⟨fk,uk h⟩Ω≤ m ∑ k=1 ∆t∥fk∥L2(Ω)∥uk h∥L2(Ω) ≤2 m ∑ k=1 ∆t∥fk∥L2(Ω)∥˜ f∥L1(L2)=2∥˜ f∥2 L1(L2)≤2∥f∥2 L1(L2). Taking the maximum in m=1,...,L, we obtain the result. Corollary 4.1. Under the same assumptions as in Lemma 4.2, we also obtain ∥˜ uh∥2 L∞(L2)+∥˜ uh∥2 L2 T(H1 0)≤C∥f∥2 L1(L2).(4.16) Now we prove some estimates for ∥∇uk h∥L3(Ω) with k=1,...,L . Corollary 4.2. For any k =1,...,L ∥∇uk h∥L3(Ω)≤C3;2h−d/6 min ∆t−1/2p2µ∥f∥L1(L2).(4.17) Proof. Let us consider the Inverse Inequality Theorem A.4 ∥∇uk h∥L3(Ω)≤C3;2h−d/6 min ∥∇uk h∥L2(Ω) ≤C3;2h−d/6 min ∆t−1/2 m ∑ k=1 ∆t∥∇uk h∥2 L2(Ω)!1/2 =C3;2h−d/6 min ∆t−1/2∥∇˜ uh∥L2(QT). Continuing this estimate with (4.15) we get the result. 100 4.2. A priori error estimates for the velocity Corollary 4.3. For any k =1,...,L ∥∇uk h∥L3(Ω)≤∆t−1/3(CShmin)−2/33 √2∥f∥2/3 L1(L2).(4.18) Proof. It is straightforward that ∥∇uk h∥L3(Ω)≤∆t−1/3∥∇˜ uh∥L3(QT) Then using the estimate (4.15) we get the result. Corollary 4.4. For any k =1,...,L ∥∇uk h∥L3(Ω)≤2C1,3;2h−1−d/6 min ∥f∥L1(L2).(4.19) Proof. We use the Inverse Inequality Theorem A.4 ∥∇uk h∥L3(Ω)≤C1,3;2h−1−d/6 min ∥uk h∥L2(Ω)≤C1,3;2h−1−d/6 min ∥˜ uh∥L∞(L2) Then using estimate (4.10) we get the result. Remark 4.5. Applying the Sobolev embedding Theorem A.2, we obtain that W1,3(Ω)⊂ L6(Ω)for d =2,3. This is, there exists a constant C6;1,3>0such that ∥u∥L6(Ω)≤C6;1,3∥∇u∥L3(Ω).(4.20) This inequality will be used along the proof of the next lemma. Lemma 4.3. Under the assumptions of Lemma 4.1 and assuming that f∈L2(QT) , it holds ∥∂∗ t˜ uh∥2 L2(QT)+1 µ∥∇˜ uh∥2 L∞(L2)+(CShmin)2∥∇˜ uh∥3 L∞(L3) ≤C1expC2h−8/3 min ∥f∥2/3 L1(L2)h−8 min∥f∥8 L2(QT)+∥f∥2 L2(QT)(4.21) where C1,C2>0depends on QTand with ∆t sufficiently small. Proof. We take the momentum equation (4.4) and we chose as test function vk h=uk h−uk−1 h: ∆t     uk h−uk−1 h ∆t     2 L2(Ω) +a(uk h,uk h−uk−1 h)+at(uk h;uk h,uk h−uk−1 h) +c(uk h,uk h,uk h−uk−1 h) = ⟨fk,uk h−uk−1 h⟩Ω. 4. Numerical approximation of transient Smagorinsky model 101 Then, ∆t     uk h−uk−1 h ∆t     2 L2(Ω) +1 µ∥∇uk h∥2 L2(Ω)+at(uk h;uk h,uk h−uk−1 h) =−c(uk h,uk h,uk h−uk−1 h)+a(uk h,uk−1 h)+⟨fk,uk h−uk−1 h⟩Ω. We sum for k=1,...,m: m ∑ k=1 ∆t     uk h−uk−1 h ∆t     2 L2(Ω) +1 µ m ∑ k=1∥∇uk h∥2 L2(Ω)+ m ∑ k=1 at(uk h;uk h,uk h−uk−1 h) =1 µ m ∑ k=1 a(uk h,uk−1 h)− m ∑ k=1 c(uk h,uk h,uk h−uk−1 h)+ m ∑ k=1⟨fk,uk h−uk−1 h⟩Ω The main idea is to compensate the right side with the left side. We will treat each term individually. • First, we use the property (A.4) to obtain 1 µ m ∑ k=1 a(uk h,uk−1 h) = m ∑ k=1 1 µ(∇uk h,∇uk−1 h)Ω ≤1 µ m ∑ k=11 2∥∇uk h∥2 L2(Ω)+1 2∥∇uk−1 h∥2 L2(Ω) =1 µ m−1 ∑ k=1∥∇uk h∥2 L2(Ω)+1 2µ∥∇um h∥2 L2(Ω). • We integrate by parts for the antisymmetric term obtaining m ∑ k=1 c(uk h,uk h,uk h−uk−1 h) = =1 2 m ∑ k=1 ∆t (uk h·∇)uk h,uk h−uk−1 h ∆t!Ω−1 2 m ∑ k=1 ∆t (uk h·∇)uk h−uk−1 h ∆t,uk h!Ω = m ∑ k=1 ∆t (uk h·∇)uk h,uk h−uk−1 h ∆t!Ω +1 2 m ∑ k=1 ∆t ∇·uk h,uk h·uk h−uk−1 h ∆t!Ω 102 4.2. A priori error estimates for the velocity We will study each term separately. Using the property (A.5) with c=1/2 we obtain that m ∑ k=1 ∆t (uk h·∇)uk h,uk h−uk−1 h ∆t!Ω ≤ m ∑ k=1 2∆t∥(uk h·∇)uk h∥2 L2(Ω)+1 8 m ∑ k=1 ∆t     uk h−uk−1 h ∆t     2 L2(Ω) . It holds ∥(uk h·∇)uk h∥2 L2(Ω)=ZΩ|(uk h·∇)uk h|dΩ2≤3ZΩ(uk h·uk h)(∇uk h:∇uk h)dΩ ≤3∥uk h∥2 L6(Ω)∥∇uk h∥2 L3(Ω) ≤3C2 6;1,3∥∇uk h∥4 L3(Ω). Henceforth, m ∑ k=1 ∆t (uk h·∇)uk h,uk h−uk−1 h ∆t!Ω ≤6C2 6;1,3 m ∑ k=1 ∆t∥∇uk h∥4 L3(Ω)+1 8 m ∑ k=1 ∆t     uk h−uk−1 h ∆t     2 L2(Ω) . On the other hand, 1 2 m ∑ k=1 ∆t ∇·uk h,uk h·uk h−uk−1 h ∆t!Ω ≤1 2 m ∑ k=1 ∆t∥∇uk h∥L3(Ω)∥uk h∥L6(Ω)     uk h−uk−1 h ∆t    L2(Ω) . Using again the same Sobolev embedding and the property (A.5) with the constant equal to 1, we obtain that 1 2 m ∑ k=1 ∆t ∇·uk h,uk h·uk h−uk−1 h ∆t!Ω ≤C6;1,3 2 m ∑ k=1 ∆t∥∇uk h∥4 L3(Ω)+1 8 m ∑ k=1 ∆t     uk h−uk−1 h ∆t     2 L2(Ω) . 4. Numerical approximation of transient Smagorinsky model 103 Summing up, m ∑ k=1 c(uk h,uk h,uk h−uk−1 h)≤(12C2 6;1,3+C6;1,3)1 2 m ∑ k=1 ∆t∥∇uk h∥4 L3(Ω) +1 4 m ∑ k=1 ∆t     uk h−uk−1 h ∆t     2 L2(Ω) . • Using again the property (A.5) with c=1, m ∑ k=1⟨fk,uk h−uk−1 h⟩Ω= m ∑ k=1 ∆t*fk,uk h−uk−1 h ∆t+Ω ≤ m ∑ k=1 ∆t∥fk∥2 L2(Ω)+1 4 m ∑ k=1 ∆t     uk h−uk−1 h ∆t     2 L2(Ω) • Using now the property (A.6), m ∑ k=1 at(uk h,uk h,uk h−uk−1 h) = m ∑ k=1∑ K∈Th (CShK)2|∇uk h|∇uk h,∇uk h−∇uk−1 hK ≥1 3∑ K∈Th (CShK)2"m ∑ k=1∥∇uk h∥3 L3(K)−∥∇uk−1 h∥3 L3(K)# =1 3∑ K∈Th (CShK)2∥∇um h∥3 L3(K)≥1 3(CShmin)2∥∇um h∥3 L3(Ω). Summing up, m ∑ k=1 ∆t     uk h−uk−1 h ∆t     2 L2(Ω) +1 µ m ∑ k=1∥∇uk h∥2 L2(Ω)+1 3(CShmin)2∥∇um h∥3 L3(Ω) ≤1 µ m−1 ∑ k=1∥∇uk h∥2 L2(Ω)+1 2µ∥∇um h∥2 L2(Ω)+C∗ 2 m ∑ k=1 ∆t∥∇uk h∥4 L3(Ω) +1 2 m ∑ k=1 ∆t     uk h−uk−1 h ∆t     2 L2(Ω) + m ∑ k=1 ∆t∥fk∥2 L2(Ω) 104 4.2. A priori error estimates for the velocity where C∗=12C2 6;1,3+C6;1,3, and we obtain that m ∑ k=1 ∆t     uk h−uk−1 h ∆t     2 L2(Ω) +1 µ∥∇um h∥2 L2(Ω)+2 3(CShmin)2∥∇um h∥3 L3(Ω) ≤C∗m ∑ k=1 ∆t∥∇uk h∥4 L3(Ω)+2 m ∑ k=1 ∆t∥fk∥2 L2(Ω)(4.22) Let us prove a bound for ∇˜ uh∈L4(0,Tf;L3(Ω)). From the previous expression, ∥∇um h∥3 L3(Ω)≤3 2(CShmin)2 C∗m ∑ k=1 ∆t∥∇uk h∥4 L3(Ω)+2 m ∑ k=1 ∆t∥fk∥2 L2(Ω)! Then, we apply the discrete Gronwall Lemma A.3 using k=∆t,B=0 and am=∥∇um h∥3 L3(Ω),bm=0,cm=2C′h−2 min∥fm∥2 L2(Ω), γm=C′C∗h−2 min∥∇um h∥L3(Ω). where C′=3 2C2 S . In order to apply this lemma, we need to guarantee that ∆tγm≤ρ with ρ<1 for all m=1,...,L. From Corollary 4.4 we obtain that ∆t∥∇um h∥L3(Ω)≤2C1,3;2∆th−1−d/6 min ∥f∥L1(L2) Then, multiplying by C′C∗h−2 min, ∆tγm≤2C1,3;2C′C∗∆th−3−d/6 min ∥f∥L1(L2)≤ρ if and only if ∆t≤h3+d/6 min 2C1,3;2C′C∗∥f∥L1(L2)·ρ for some ρ∈(0,1) . We select Corollary 4.4 since it seems to be the least restrictive compared to Corollary 4.2 and 4.3. Now, applying the discrete Gronwall Lemma A.3, ∥∇um h∥L3(Ω)≤exp m ∑ k=1 ∆tγk 1−∆tγk!· 2C′h−2 min m ∑ k=1 ∆t∥fk∥2 L2(Ω)! ≤exp 1 1−ρ m ∑ k=1 ∆tγk!· 2C′h−2 min m ∑ k=1 ∆t∥fk∥2 L2(Ω)! =exp C′C∗h−2 min 1−ρ m ∑ k=1 ∆t∥∇uk h∥L3(Ω)!· 2C′h−2 min m ∑ k=1 ∆t∥fk∥2 L2(Ω)! 4. Numerical approximation of transient Smagorinsky model 105 And taking the supreme in m=1,...,L, ∥∇˜ uh∥L∞(L3)≤exp C′C∗h−2 min 1−ρ∥∇˜ uh∥L1(L3)!2C′h−2 min∥f∥2 L2(QT) The exponential term is bounded since ∇˜ uh∈L3(0,Tf;L3(Ω)) from (4.15) and L3(0,Tf),→L1(0,Tf), then, there exists C′ Tf>0 such that ∥∇˜ uh∥L1(L3)≤C′ Tf∥∇˜ uh∥L3(QT)≤ 3 √2C′ Tf (CShmin)2/3∥f∥2/3 L1(L2) with this, we have proven that ˜ uh∈L∞(0,Tf,W1,3(Ω)). Again, since L∞(0,Tf),→L4(0,Tf), there exists C⋆ Tf>0 such that ∥∇˜ uh∥L4(L3)≤C⋆ Tf∥∇˜ uh∥L∞(L3). Then, taking the supreme in m=1,...,Lin (4.22) we obtain ∥∂∗ t˜ uh∥2 L2(QT)+1 µ∥∇˜ uh∥2 L∞(L2)+(CShmin)22 3∥∇˜ uh∥3 L∞(L3) ≤C∗∥∇˜ uh∥4 L4(L3)+2∥f∥2 L2(QT) ≤C∗(C∗ Tf)4∥∇˜ uh∥4 L∞(L3)+∥f∥2 L2(QT) ≤C∗(C∗ Tf)4exp 4C′C∗h−2 min 1−ρ∥∇˜ uh∥L1(L3)!2C′h−2 min∥f∥2 L2(QT)4+∥f∥2 L2(QT) ≤16C∗(C⋆ TfC′)4exp  43 √2C′C∗C′ Tfh−8/3 min C2/3 S(1−ρ)∥f∥2/3 L1(L2) h−8 min∥f∥8 L2(QT)+∥f∥2 L2(QT) and we obtain the result. Corollary 4.5. Under the same assumptions as in Lemma 4.3, we also obtain ∥∂∗ t˜ uh∥2 L2(QT)+∥˜ uh∥2 L∞ T(H1 0) ≤C1expC2h−8/3 min ∥f∥2/3 L1(L2)h−8 min∥f∥8 L2(QT)+∥f∥2 L2(QT)(4.23) where C1,C2>0depends on QTand ∆t sufficiently small. Remark 4.6. We are working with the Smagorinsky model, therefore we have set a mesh Th with h>0 . The Lemma 4.3 shows an estimate with negative power on hmin which should not be a issue since hmin =0. 112 5.1. Problem statement where the eddy viscosity term is given by νt(˜ w) = ∑ K∈Th (CShK)2|∇˜ w|K|XK. To simplify notation, we introduce the space ˜ Xh=˜ Yhט Yhט Qh , this is, the space formed by piecewise constant functions ˜ Vh related to the velocity field ˜ vh∈˜ Yh , its discrete time derivative ∂∗ t˜ vh∈˜ Yh and pressure ˜qh∈˜ Qh . In order to consider nonhomogeneous Dirichlet BCs for the application of the RB method, we set up a divergence free lift uD that will take part of the formulation, then, ˜ wh=˜ uh+uD where ˜ uh verifies homogeneous Dirichlet BCs. We assume that uD∈L∞(0,Tf;H1 0(Ω)) and for simplicity, it must verify that uD(0,·) is equal to the initial condition, then, ˜ wh(0,·) = 0 . We consider µ=Re the parameter of the problem and D the parameter set. Then the problem to solve in its compact form is defined as follows (Given µ∈Dfind ˜ Uh(µ) = (˜ uh(µ),˜ph(µ)) ∈˜ Xhsuch that ˜ A(˜ Uh(µ),˜ Vh;µ) = F(˜ Vh;µ),∀˜ Vh∈˜ Xh,(5.1) where A is defined as ˜ A(˜ Uh(µ),˜ Vh;µ) = ˜ A0(˜ Uh(µ),˜ Vh;µ)+ ˜ A1(˜ Uh(µ),˜ Vh;µ)+ ˜ A2(˜ Uh(µ),˜ Vh;µ)(5.2) with ˜ A0(˜ U,˜ V;µ) = ZQT ∂∗ t˜ u·˜ vdQT ˜ A1(˜ U,˜ V;µ) = 1 µZQT ∇˜ u:∇˜ vdQT+1 2ZQT [(˜ u·∇)˜ u·˜ v−(˜ u·∇)˜ v·˜ u]dQT +1 2ZQT [(uD·∇)˜ u·˜ v−(uD·∇)˜ v·˜ u]dQT +1 2ZQT [(˜ u·∇)uD·˜ v−(˜ u·∇)˜ v·uD]dQT +ZQT νt(˜ u+uD)∇(˜ u+uD):∇˜ vdQT ˜ A2(˜ U,˜ V;µ) = ZQT (∇·˜ u)˜qdQT−ZQT (∇·˜ v)˜pdQT F(˜ V;µ) = ZQT f·˜ vdQT−1 µZQT ∇uD:∇˜ vdQT −1 2ZQT (uD·∇)uD·˜ vdQT+1 2ZQT (uD·∇)˜ v·uDdQT. (5.3) 5. A posteriori error estimation for unsteady Smagorinsky model 113 where ˜ U= (˜ u,∂∗ t˜ u,˜p) and ˜ V= (˜ v,∂∗ t˜ v,˜q) . Note that ˜ A0 is the operator related to the derivative on time, ˜ A1 involves convection and viscous effects and ˜ A2 involves the pressure. From now on, we consider the following norm on ˜ Xh : ∥˜ Vh∥˜ X=q∥˜ vh∥2 L2 T(H1 0)+∥˜qh∥2 L2(QT),(5.4) where we recall that the Lp T(H1 0(Ω)) -norm was introduced in (4.13) - (4.14) for p∈[1,+∞] . With this setting, we can state the Reduced Basis problem as (Given µ∈D, find ˜ UN(µ)∈˜ XNsuch that ˜ A(˜ UN(µ),˜ VN;µ) = F(˜ VN;µ),∀˜ VN∈˜ XN (5.5) where ˜ XN=˜ YNט YNט QN , with ˜ YN⊂˜ Yh and ˜ QN⊂˜ Qh are the reduced basis spaces. Remark 5.1. The results of existence and stability proved in Chapter 4 still hold for problem (5.1) for smooth enough lift uD. 5.2 Well-posedness analysis To study the well-posedness of the model, we use the Brezzi-Rappaz-Raviart (BRR) theory (see, e.g. [ 8 ]). This theory provides a technique to obtain an a posteriori error estimator for nonlinear parametric problems and it has already applied to the steady case. We study the directional derivative defined in Definition A.1 of operator ˜ A to apply this theory. Let ˜ Z= (˜ z,∂∗ t˜ z,˜r) , it holds, ∂1˜ A0(˜ U,˜ V;µ)( ˜ Z) =ZQT ∂∗ t˜ z·˜ vdQT, ∂1˜ A1(˜ U,˜ V;µ)( ˜ Z) = 1 µZQT ∇˜ z:∇˜ vdQT +1 2ZQT [( ˜ w·∇)˜ z·˜ v−(˜ w·∇)˜ v·˜ z]dQT +1 2ZQT [(˜ z·∇)˜ w·˜ v−(˜ z·∇)˜ v·˜ w]dQT +ZQT νt(˜ w)∇˜ z:∇˜ vdQT +∑ K∈ThZTf 0ZK (CShK)2∇(˜ w):∇˜ z |∇(˜ w)|∇(˜ w):∇˜ vdKdt, ∂1˜ A2(˜ U,˜ V;µ)( ˜ Z) =ZQT (∇·˜ z)˜qdQT−ZQT (∇·˜ v)˜rdQT (5.6) where ˜ w=˜ u+uD . 114 5.2. Well-posedness analysis Remark 5.2. Before to start, we define some constants derived from relevant results presented in Appendix A.2: • The norm related to the inner product (2.12) is equivalent to the H1 0(Ω) -norm, then, there exists a constant C1,2;T>0such that ∥∇˜ vh∥L2(Ω)≤C1,2;T∥˜ vh∥T,∀˜ vh∈˜ Yh,∀t∈If.(5.7) • Applying the Sobolev embedding Theorem A.2 we obtain that H1 0(Ω),→L4(Ω),→ L2(Ω) with constants C4;1,2,C2;1,2>0 . Thanks to the previous point, there exist constants C4;T=C4;1,2C1,2;Tand C2;T=C2;1,2C1,2;Tsuch that ∥˜ vh∥L4(Ω)≤C4;T∥v∥T,∥˜ vh∥L2(Ω)≤C2;T∥v∥T,∀˜ vh∈˜ Yh,∀t∈If.(5.8) • For all K∈Th , we apply the Local Inverse Inequality Theorem A.3 as in (A.11) , introducing the constant C3;2 >0, this is, ∥∇˜ vh∥L3(K)≤C3;2h−d/6 K∥∇˜ vh∥L2(K),∀K∈Th,∀˜ vh∈˜ Yh,∀t∈If.(5.9) We start proving the continuity of operator ∂1˜ A . Proposition 5.1. For any µ∈Dand ˜ Uh(µ)∈˜ Xh, it holds |∂1˜ A(˜ Uh(µ),˜ Vh;µ)( ˜ Zh)|≤γh∥˜ Zh∥˜ X∥˜ Vh∥˜ X,∀˜ Zh,˜ Vh∈˜ Xh,(5.10) where γh=maxµ∈Dγh(µ)and γh(µ) = 2C2 2;T ∆t+2C1,2;T+C2 1,2;T µ+2C1,2;TC2 4;T∥˜ wh∥L∞ T(H1 0) +2(C1,2;TC3;2CS)2h2−d/3 max ∥∇˜ wh∥L∞(L3). with ˜ wh=˜ uh+uD , and the constants C1,2;T , C2;T , C4;T and C3;2 are defined in (5.7) - (5.9) . Proof. We start from |∂1˜ A(˜ Uh(µ),˜ Vh;µ)( ˜ Zh)|≤|∂1˜ A0(˜ Uh(µ),˜ Vh;µ)( ˜ Zh)| +|∂1˜ A1(˜ Uh(µ),˜ Vh;µ)( ˜ Zh)|+|∂1˜ A2(˜ Uh(µ),˜ Vh;µ)( ˜ Zh)| We study each term separately using the constants defined in (5.7)-(5.9). 5. A posteriori error estimation for unsteady Smagorinsky model 115 The first term involves the derivative in time, |∂1˜ A0(˜ Uh(µ),˜ Vh)( ˜ Zh)|= L ∑ k=1 ∆tZΩ zk h−zk−1 h ∆t·vk hdΩ ≤ L ∑ k=1 ∆t     zk h−zk−1 h ∆t    L2(Ω)∥vk h∥L2(Ω)≤C2 2;T ∆t L ∑ k=1 ∆t∥zk h−zk−1 h∥T∥vk h∥T ≤2C2 2;T ∆t∥˜ zh∥L2 T(H1 0)∥˜ vh∥L2 T(H1 0)≤2C2 2;T ∆t∥˜ Zh∥˜ X∥˜ Vh∥˜ X (5.11) Next, we study |∂1˜ A1(˜ Uh(µ),˜ Vh;µ)( ˜ Zh)|term by term,  1 µZQT ∇˜ zh:∇˜ vhdQT≤1 µ L ∑ k=1 ∆tZΩ ∇zk h:∇vk hdΩ ≤1 µ L ∑ k=1 ∆t∥∇zk h∥L2(Ω)∥∇vk h∥L2(Ω)≤C2 1,2;T µ L ∑ k=1 ∆t∥zk h∥T∥vk h∥T Furthermore,  1 2ZQT [( ˜ wh·∇)˜ zh·˜ vh−(˜ wh·∇)˜ vh·˜ zh]dQT +1 2ZQT [(˜ zh·∇)˜ wh·˜ vh−(˜ zh·∇)˜ vh·˜ wh]dQT ≤1 2 L ∑ k=1 ∆tZΩ(wk h·∇)zk h·vk h−(wk h·∇)vk h·zk hdΩ +1 2 L ∑ k=1 ∆tZΩ(zk h·∇)wk h·vk h−(zk h·∇)vk h·wk hdΩ ≤1 2 L ∑ k=1 ∆th∥wk h∥L4(Ω)∥∇zk h∥L2(Ω)∥vk h∥L4(Ω)+∥wk h∥L4(Ω)∥∇vk h∥L2(Ω)∥zk h∥L4(Ω) +∥zk h∥L4(Ω)∥∇wk h∥L2(Ω)∥vk h∥L4(Ω)+∥zk h∥L4(Ω)∥∇vk h∥L2(Ω)∥wk h∥L4(Ω)i ≤2C1,2;TC2 4;T L ∑ k=1 ∆t∥wk h∥T∥zk h∥T∥vk h∥T. where we use (5.7) and (5.8) and the Hölder’s Inequality A.1. 116 5.2. Well-posedness analysis Finally, the viscosity term is bounded by ZQT νt(˜ wh)∇˜ zh:∇˜ vhdQT+∑ K∈ThZTf 0∆tZK (CShK)2∇˜ wh:∇˜ zh |∇˜ wh|∇˜ wh:∇˜ vhdKdt ≤2∑ K∈Th L ∑ k=1 ∆tZK (CShK)2∥∇wk h∥L3(K)∥∇zk h∥L3(K)∥∇vk h∥L3(K)dK ≤2C2 3;2 ∑ K∈Th L ∑ k=1 ∆tZK C2 Sh2−d/3 K∥∇wk h∥L3(K)∥∇zk h∥L2(K)∥∇vk h∥L2(K)dK ≤2(C3;2CS)2h2−d/3 max L ∑ k=1 ∆th∥∇wk h∥L3(Ω)∥∇zk h∥L2(Ω)∥∇vk h∥L2(Ω)i ≤2(C1,2;TC3;2CS)2h2−d/3 max L ∑ k=1 ∆th∥∇wk h∥L3(Ω)∥zk h∥T∥vk h∥Ti where ˜ wh=˜ uh+uD and the constants C1,2;T is defined in (5.7) and C3;2 comes from the application of the Local Inverse Inequality (5.9). To sum up, |∂1˜ A1(˜ Uh(µ),˜ Vh;µ)( ˜ Zh)| ≤ L ∑ k=1 ∆t C2 1,2;T µ+2C1,2;TC2 4;T∥wk h∥T +2(C1,2;TC3;2CS)2h2−d/3 max ∥∇wk h∥L3(Ω)∥zk h∥T∥vk h∥T ≤ C2 1,2;T µ+2C1,2;TC2 4;T∥˜ wh∥L∞ T(H1 0) +2(C1,2;TC3;2CS)2h2−d/3 max ∥∇˜ wh∥L∞(L3)∥˜ zh∥L2 T(H1 0)∥˜ vh∥L2 T(H1 0) ≤ C2 1,2;T µ+2C1,2;TC2 4;T∥˜ wh∥L∞ T(H1 0) +2(C1,2;TC3;2CS)2h2−d/3 max ∥∇˜ wh∥L∞(L3)∥˜ Zh∥˜ X∥˜ Vh∥˜ X (5.12) Finally, |∂1˜ A2(˜ Uh(µ),˜ Vh)( ˜ Zh)|=ZQT (∇·˜ zh)˜qhdQT−ZQT (∇·˜ vh)˜rhdQT ≤ZTf 0∥∇˜ zh∥L2(Ω)∥˜qh∥L2(Ω)dt +ZTf 0∥∇˜ vh∥L2(Ω)∥˜rh∥L2(Ω)dt ≤C1,2;T∥˜ zh∥L2 T(H1 0)∥˜qh∥L2(QT)+C1,2;T∥˜ vh∥L2 T(H1 0)∥˜rh∥L2(QT) ≤2C1,2;T∥˜ Zh∥˜ X∥˜ Vh∥˜ X 5. A posteriori error estimation for unsteady Smagorinsky model 117 We get the result with γh(µ) = 2C2 2;T ∆t+2C1,2;T+C2 1,2;T µ+2C1,2;TC2 4;T∥˜ wh∥L∞ T(H1 0) +2(C1,2;TC3;2CS)2h2−d/3 max ∥∇˜ wh∥L∞(L3). and we set γh=max µ∈Dγh(µ). The continuity constant γh grows as ∆t−1 because of the derivative in time. This dependency is admissible for moderate time steps. If we include the time derivative in the ˜ X -norm, the constant continuity does not depend on the time step. However, the dependency on ∆t−1 appears when proving the inf-sup condition. Proposition 5.2. For any µ∈D and ˜ Uh(µ)∈˜ Xh , let us suppose that ∥uD∥L∞ T(H1 0)<1 C1,2;TC2 4;T and ∥˜ uh∥L∞ T(H1 0)<1 C1,2;TC2 4;T−∥uD∥L∞ T(H1 0). where the constants C1,2;T and C4;T are defined in (5.7) and (5.8) Then, for all µ∈D it holds, βh(µ) = inf ˜ Zh∈˜ Xh sup ˜ Vh∈˜ Xh ∂1˜ A(˜ Uh(µ),˜ Vh;µ)( ˜ Zh) ∥˜ Zh∥˜ X∥˜ Vh∥˜ X ,(5.13) where βh(µ) = Ct+(1+γ2 h) 1+γ2 h ˜ β2 h(µ)+2√CCt ˜ βh(µ)!!−1/2 with ˜ βh(µ) = 1−C1,2;TC2 4;T∥˜ wh∥L∞ T(H1 0),Ct=Tf ∆t and γhis the continuity constant deduced in Proposition 5.1 that grows as ∆t−1. Proof. We base this proof on Theorem 1 in [ 10 ]. The key is to prove that for a given µ∈D, and ˜ Zh∈˜ Xh, ∥˜ Zh∥˜ Xβh(µ)≤˜ S(˜ Uh(µ),˜ Zh) = sup ˜ Vh∈˜ Xh ∂1˜ A(˜ Uh(µ),˜ Vh;µ)( ˜ Zh) ∥˜ Vh∥˜ X . 118 5.2. Well-posedness analysis First, we start deducing an estimate for the velocity in its related norm. We select ˜ Vh=˜ Zh∈Xh, it holds ∂1˜ A(˜ Uh(µ),˜ Zh;µ)( ˜ Zh) = ∂1˜ A0(˜ Uh(µ),˜ Zh)( ˜ Zh) +∂1˜ A1(˜ Uh(µ),˜ Zh;µ)( ˜ Zh)+∂1˜ A2(˜ Uh(µ),˜ Zh)( ˜ Zh). The time derivative term is positive since ∂1˜ A0(˜ Uh(µ),˜ Zh)( ˜ Zh) = ZQT ∂∗ t˜ zh·˜ zhdQT= L ∑ k=1 ∆t zk h−zk−1 h ∆t,zk h! =1 2 L ∑ k=1∥zk h∥2 L2(Ω)−∥zk−1 h∥2 L2(Ω)+∥zk h−zk−1 h∥2 L2(Ω) =1 2∥zL h∥2 L2(Ω)+1 2 L ∑ k=1∥zk h−zk−1 h∥2 L2(Ω); For the second term, on one side, we bound the terms related to viscosity, ZQT1 µ+νt(˜ wh)|∇˜ zh|2dQT+∑ K∈ThZTf 0ZK (CShK)2|∇˜ wh:∇˜ zh|2 |∇˜ wh|dΩdt ≥ZQT1 µ+νt(˜ wh)|∇˜ zh|2dQT≥∥˜ zh∥2 L2 T(H1 0); where we recall that ˜ wh=˜ uh+uD. On another hand, 1 2ZQT [(˜ zh·∇)˜ wh·˜ zh−(˜ zh·∇)˜ zh·˜ wh]dQT =1 2 L ∑ k=1 ∆tZΩ[(zk h·∇)wk h·zk h−(zk h·∇)zk h·wk h]dΩ ≤1 2 L ∑ k=1 ∆t[∥zk h∥L4(Ω)∥∇wk h∥L2(Ω)∥zk h∥L4(Ω)+∥zk h∥L4(Ω)∥∇zk h∥L2(Ω)∥wk h∥L4(Ω)] ≤C1,2;TC2 4;T L ∑ k=1 ∆t∥wk h∥T∥zk h∥2 T≤C1,2;TC2 4;T∥˜ wh∥L∞ T(H1 0)∥˜ zh∥2 L2 T(H1 0). Lastly, ∂1˜ A2(˜ Uh(µ),˜ Zh)( ˜ Zh) = 0, then, ∂1˜ A(˜ Uh(µ),˜ Zh;µ)( ˜ Zh)≥(1−C1,2;TC2 4;T∥˜ wh∥L∞ T(H1 0))∥˜ zh∥2 L2 T(H1 0).(5.14) 5. A posteriori error estimation for unsteady Smagorinsky model 119 We need that 1−C1,2;TC2 4;T∥˜ wh∥L∞ T(H1 0)>0 and that is equivalent to 1 C1,2;TC2 4;T >∥˜ uh∥L∞ T(H1 0)+∥uD∥L∞ T(H1 0) Then, if we suppose that ∥uD∥T<1 C1,2;TC2 4;T , we need that ∥˜ uh∥L∞ T(H1 0)<1 C1,2;TC2 4;T−∥uD∥L∞ T(H1 0). Note that inequality (5.14) is possible thanks to Corollary 4.5 for small data. Finally, there exists ˜ βh(µ) = 1−C1,2;TC2 4;T∥˜ wh∥L∞ T(H1 0)>0 such that ∥˜ zh∥2 L2 T(H1 0)˜ βh(µ)≤˜ S(˜ Uh(µ),˜ Zh)∥˜ Zh∥˜ X,∀˜ Zh∈˜ Xh. Next, we deduce the pressure estimate. Setting Uk h,Vk h,Zk h∈˜ Xh with Xh=Yh×Yh×Qh defined by Uk h= (uk h,uk−1 h,pk h),Vk h= (vk h,vk−1 h,qk h),Zk h= (zk h,zk−1 h,rk h), from definition (5.6), we can express ∂1˜ A(˜ Uh(µ),˜ Vh;µ)( ˜ Zh) = L ∑ k=1 ∆t∂1A(Uk h(µ),Vk h;µ)(Zk h)(5.15) where ∂1A(Uk h(µ),Vk h;µ)(Zk h) = ∂1A0(Uk h(µ),Vk h;µ)(Zk h)+∂1A1(Uk h(µ),Vk h;µ)(Zk h) +∂1A2(Uk h(µ),Vk h;µ)(Zk h). with the forms ∂1A1and ∂1A2as the same as ∂1˜ A1and ∂1˜ A2by changing the integration domain to Ωinstead of QTand ∂1A0(Uk h(µ),Vk h;µ)(Zk h) = ZΩ zk h−zk−1 h ∆t·vk hdΩ. Choosing the test function Vk h= (vk h,vk−1 h,0), for all k=1,...,L, we have ∂1A(Uk h(µ),Vk h;µ)(Zk h) = ∂1A0(Uk h(µ),Vk h;µ)(Zk h) +∂1A1(Uk h(µ),Vk h;µ)(Zk h)−ZΩ(∇·vk h)rk hdΩ 120 5.2. Well-posedness analysis from where ZΩ(∇·vk h)rk hdΩ=−∂1A(Uk h(µ),Vk h;µ)(Zk h) +∂1A0(Uk h(µ),Vk h;µ)(Zk h)+∂1A1(Uk h(µ),Vk h;µ)(Zk h) Now, we define a supreme Skrelated to each time step, for a given µ∈D: Sk(Uk h(µ),Zk h) = sup Vk h∈Xh ∂1A(Uk h(µ),Vk h;µ)(Zk h) ∥Vk h∥X ,∀Zk h∈Xh,(5.16) with ∥Vk h∥X=q∥vk h∥2 T+∥qk h∥2 L2(Ω),∀Vk h∈Xh. Because of the Vk hselection, we obtain that ∂1A(Uk h(µ),Vk h;µ)(Zk h)≤Sk(Uk h(µ),Zk h)∥vk h∥T,∀Zk h∈Xh Analogously to (5.11) and (5.12), it can be obtained: ∂1A0(Uk h(µ),Vk h;µ)(Zk h)≤C2 2;T ∆t∥zk h−zk−1 h∥T∥vk h∥T, ∂2A0(Uk h(µ),Vk h;µ)(Zk h)≤γk h(µ)∥zk h∥T∥vk h∥T, where γk h(µ) = C2 1,2;T µ+2C1,2;TC2 4;T∥wk h∥T+2(C1,2;TC3;2CS)2h2−d/3 max ∥∇wk h∥L3(Ω). with wk h=uk h+uD. Then, ZΩ(∇·vk h)rk hdΩ≤Sk(Uk h(µ),Zk h)∥vk h∥T+C2 2;T ∆t∥zk h−zk−1 h∥T∥vk h∥T+γk h(µ)∥zk h∥T∥vk h∥T Dividing by ∥vk h∥T , we are able to apply the Inf-sup condition Theorem A.8 since the T-norm is equivalent to the H1 0(Ω)-norm. There exists α>0 such that α∥rk h∥L2(Ω)≤Sk(Uk h(µ),Zk h)+ C2 2;T ∆t∥zk h−zk−1 h∥T+γ∗ h∥zk h∥T. where γ∗ h=max µ∈D,k=1,...,Lγk h(µ). Then, α2∥rk h∥2 L2(Ω)≤3 S2 k(Uk h(µ),Zk h)+ C4 2;T ∆t2∥zk h−zk−1 h∥2 T+(γ∗ h)2∥zk h∥2 T!, 5. A posteriori error estimation for unsteady Smagorinsky model 121 and multiplying by ∆tand summing on k, α2∥˜rh∥2 L2(QT)≤3 L ∑ k=1 ∆t"S2 k(Uk h(µ),Zk h)+3C4 2;T ∆t2∥zk h−zk−1 h∥2 T+3(γ∗ h)2∥zk h∥2 T#, We need an upper bound of this inequality in terms of the supreme ˜ S2(˜ Uh(µ),˜ Zh) and ∥˜ zh∥2 L2 T(H1 0). On one side, α2∥˜rh∥2 L2(QT)≤3 L ∑ k=1 ∆thS2 k(Uk h(µ),Zk h)i+3 2C4 2;T ∆t2+(γ∗ h)2!∥˜ zh∥2 L2 T(H1 0) ≤3 L ∑ k=1 ∆thS2 k(Uk h(µ),Zk h)i+6γ2 h∥˜ zh∥2 L2 T(H1 0), where γh is the continuity constant deduced in Proposition 5.1. On the other side, let us define the functions ˜ Vk h∈˜ Xhfor all k=1,...,Lby ˜ Vk h=(Vk hif t∈(tk−1,tk), 0 otherwise, Then, using the definition of ˜ S(˜ Uh(µ),˜ Zh), ˜ S(˜ Uh(µ),˜ Zh) = sup ˜ Vh∈˜ Xh ∂1˜ A(˜ Uh(µ),˜ Vh;µ)( ˜ Zh) ∥˜ Vh∥˜ X≥sup ˜ Vk h∈˜ Xh ∂1˜ A(˜ Uh(µ),˜ Vk h;µ)( ˜ Zh) ∥˜ Vk h∥˜ X , for all k=1,...,L. Note that ˜ Vk hbelongs to ˜ Xhand ∥˜ Vk h∥2 ˜ X= L ∑ n=1 ∆t∥vn h∥2 T+∥qn h∥2 L2(Ω)=∆t∥vk h∥2 T+∥qk h∥2 L2(Ω)=∆t∥Vk h∥2 X with Vk h∈Xh. Moreover, note that from (5.15), ∂1˜ A(˜ Uh(µ),˜ Vk h;µ)( ˜ Zh) = ∆t∂1A(Uk h(µ),Vk h;µ)(Zk h) then for all k=1,...,L, ˜ S(˜ Uh(µ),˜ Zh)≥sup Vk h∈Xh ∆t∂1A(Uk h(µ),Vk h;µ)(Zk h) √∆t∥Vk h∥X =√∆tSk(Uk h(µ),Zk h). in particular, L ∑ k=1 ∆tS2 k(Uk h(µ),Zk h)≤ L ∑ k=1 ˜ S(˜ Uh(µ),˜ Zh) = L˜ S(˜ Uh(µ),˜ Zh) = Tf ∆t ˜ S(˜ Uh(µ),˜ Zh) 128 5.3. A posteriori error estimator Then, thanks to (5.18), we can write ⟨DA(˜ UN(µ);µ)(H(˜ Z1 h;µ)−H(˜ Z2 h;µ)),˜ Vh⟩≤ρT∥˜ UN(µ)−˜ ξ∥˜ X∥˜ Z1 h−˜ Z2 h∥˜ X∥˜ Vh∥˜ X (5.28) From (5.22) ⟨DA(˜ UN(µ);µ)(H(˜ Z1 h;µ)−H(˜ Z2 h;µ)),Tµ N(H(˜ Z1 h;µ)−H(˜ Z2 h;µ))⟩ ≥βN(µ)∥H(˜ Z1 h;µ)−H(˜ Z2 h;µ)∥˜ X∥Tµ N(H(˜ Z1 h;µ)−H(˜ Z2 h;µ))∥˜ X And applying (5.28) to ˜ Vh=Tµ N(H(˜ Z1 h;µ)−H(˜ Z2 h;µ)) we have that, βN(µ)∥H(˜ Z1 h;µ)−H(˜ Z2 h;µ)∥˜ X∥Tµ N(H(˜ Z1 h;µ)−H(˜ Z2 h;µ))∥˜ X ≤ρT∥˜ UN(µ)−˜ ξ∥˜ X∥˜ Z1 h−˜ Z2 h∥˜ X∥Tµ N(H(˜ Z1 h;µ)−H(˜ Z2 h;µ))∥˜ X. Note that Tµ N˜ Zh=0 if and only if ˜ Zh=0 . Therefore, Tµ N(H(˜ Z1 h;µ)−H(˜ Z2 h;µ)) =0 if H(˜ Z1 h;µ)=H(˜ Z2 h;µ). Simplifying ∥H(˜ Z1 h;µ)−H(˜ Z2 h;µ)∥˜ X≤ρT βN(µ)∥˜ UN(µ)−˜ ξ∥˜ X∥˜ Z1 h−˜ Z2 h∥˜ X. If ˜ Z1 hand ˜ Z2 hare in B˜ X(˜ UN(µ),α), then ∥˜ UN(µ)−˜ ξ∥˜ X≤αand ∥H(˜ Z1 h;µ)−H(˜ Z2 h;µ)∥˜ X≤ρT βN(µ)α∥˜ Z1 h−˜ Z2 h∥˜ X. Then, H(·;µ) is a contraction if α<βN(µ) ρT . So it follows that there can exist at most one fixed point of H(·;µ) inside B˜ X˜ UN(µ),βN(µ) ρT , and hence, at most one solution ˜ Uh(µ)of (5.1) in this ball. Now we are in a position to define an a posteriori error bound estimator by ∆N(µ) = βN(µ) 2ρTh1−p1−τN(µ)i,∀µ∈D,(5.29) where βN(µ) is defined in (5.19) , ρT is defined in (5.18) and τN(µ) = 4εN(µ)ρT βN(µ)2,(5.30) with εN(µ) = ∥R(˜ UN(µ);µ)∥˜ X′=sup ˜ Vh∈˜ Xh ⟨R(˜ UN(µ);µ),˜ Vh⟩ ∥˜ Vh∥˜ X .(5.31) We observe that τN(µ) in (5.30) is a re-scaling of the residual (5.31) . Note that ∆N(µ) in (5.29) is only defined if τN(µ)≤1 , which holds only if the residual is sufficiently small. The error estimator is given by the following theorem. 5. A posteriori error estimation for unsteady Smagorinsky model 129 Theorem 5.2 (Existence).Assume that βN(µ)>0 and τN(µ)≤1 for all µ∈D . Then there exists a unique solution ˜ Uh(µ) of (5.1) such that the error with respect ˜ UN(µ) , solution of (5.5), is bounded by the a posteriori error bound estimator, i.e., ∥˜ Uh(µ)−˜ UN(µ)∥˜ X≤∆N(µ),(5.32) with effectivity, ∆N(µ)≤2γN(µ) βN(µ)+τN(µ)∥˜ Uh(µ)−˜ UN(µ)∥˜ X.(5.33) Proof. The main goal of this proof is to show that there exists a unique ˜ Uh(µ)∈˜ Xh close to ˜ UN(µ)∈˜ XNsolution of (5.1) and (5.5) (resp.), i.e., ⟨R(˜ Uh(µ);µ),˜ Vh⟩=0,∀˜ Vh∈˜ Xh.(5.34) To do this, it is enough to prove that H is a contraction in a certain compact subset K of ˜ Xh such that H(K)⊂K . Therefore, H has a fixed point thanks to the fixed-point theorem [ 9 ]. The operator H is continuous by the Proposition 5.1, then we only need to prove that it is contractive. We consider the definition at the beginning of the previous proof. Then, H(˜ Zh;µ)−˜ UN(µ) = ˜ Zh−˜ UN(µ)−DA(˜ UN(µ);µ)−1R(˜ Zh;µ). Multiplying by DA(˜ UN(µ);µ), we obtain ⟨DA(˜ UN(µ);µ)(H(˜ Zh;µ)−˜ UN(µ)),˜ Vh⟩=⟨DA(˜ UN(µ);µ)( ˜ Zh−˜ UN(µ)),˜ Vh⟩ −⟨R(˜ Zh;µ)−R(˜ UN(µ);µ),˜ Vh⟩−⟨R(˜ UN(µ);µ),˜ Vh⟩,∀˜ Vh∈˜ Xh As in the proof of Theorem 5.1, it holds R(˜ Zh;µ)−R(˜ UN(µ);µ) = DA(˜ ξ(µ);µ)( ˜ Zh−˜ UN(µ)), where ˜ ξ(µ) = s˜ Zh+(1−s)˜ UN(µ) , s∈(0,1) . By this way, and thanks to Lemma 5.1, we obtain: ⟨DA(˜ UN(µ);µ)(H(˜ Zh;µ)−˜ UN(µ)),˜ Vh⟩ =⟨DA(˜ UN(µ);µ)( ˜ Zh−˜ UN(µ)),˜ Vh⟩ −⟨DA(˜ ξ(µ);µ)( ˜ Zh−˜ UN(µ)),˜ Vh⟩−⟨R(˜ UN(µ);µ),˜ Vh⟩ =⟨(DA(˜ UN(µ);µ)−DA(˜ ξ(µ);µ))( ˜ Zh−˜ UN(µ)),˜ Vh⟩ −⟨R(˜ UN(µ);µ),˜ Vh⟩ ≤ρT∥˜ UN(µ)−˜ ξ(µ)∥˜ X∥˜ Zh−˜ UN(µ)∥˜ X∥˜ Vh∥˜ X+εN(µ)∥˜ Vh∥˜ X ≤(ρT∥˜ Zh−˜ UN(µ)∥2 ˜ X+εN(µ))∥˜ Vh∥˜ X. 130 5.3. A posteriori error estimator Then, using (5.22) βN(µ)∥H(˜ Zh;µ)−˜ UN(µ)∥˜ X∥Tµ N(H(˜ Zh;µ)−˜ UN(µ))∥˜ X ≤⟨DA(˜ UN(µ);µ)(H(˜ Zh;µ)−˜ UN(µ)),Tµ N(H(˜ Zh;µ)−˜ UN(µ))⟩ ≤(ρT∥˜ Zh−˜ UN(µ)∥2 ˜ X+εN(µ))∥Tµ N(H(˜ Zh;µ)−˜ UN(µ))∥˜ X Then as ˜ Zh∈B˜ X(˜ UN(µ),α)and simplifying, ∥H(˜ Zh;µ)−˜ UN(µ)∥˜ X≤ρT βN(µ)α2+εN(µ) βN(µ) In order to ensure that H maps B˜ X(˜ UN(µ),α) into a part of itself, we are seeking the values of αsuch that ρT βN(µ)α2+εN(µ) βN(µ)≤α This holds if αis between the two roots of the second order equation ρTα2−βN(µ)α+εN(µ) = 0 which are, α±=βN(µ)±pβN(µ)2−4ρTεN(µ) 2ρT =βN(µ) 2ρT"1±s1−4ρTεN(µ) βN(µ)2#. Observe that as τN(µ)≤1 , then α−≤α+≤βN(µ) ρT . Consequently, if α−≤α≤α+ , there exists a unique solution ˜ Uh(µ) to (5.1) in the ball B˜ X(˜ UN(µ),α) . To obtain the estimator, we take α=α−=∆N(µ). To prove the efficiency, let us define the error ˜ Eh(µ) = ˜ Uh(µ)−˜ UN(µ) . From the definition of the residual and applying the mean value theorem, for some s∈(0,1) we have that ⟨R(˜ UN(µ);µ),˜ Vh⟩=˜ A(˜ UN(µ),˜ Vh;µ)−F(˜ Vh;µ) =˜ A(˜ UN(µ),˜ Vh;µ)−˜ A(˜ Uh(µ),˜ Vh;µ) =∂1˜ A(s˜ Uh(µ)+(1−s)˜ UN(µ),˜ Vh;µ)( ˜ Eh(µ)) =⟨DA(s˜ Uh(µ)+(1−s)˜ UN(µ);µ)˜ Eh(µ),˜ Vh⟩ =⟨(DA(s˜ Uh(µ)+(1−s)˜ UN(µ);µ)−DA(˜ UN(µ);µ)) ˜ Eh(µ),˜ Vh⟩ +⟨DA(˜ UN(µ);µ)˜ Eh(µ),˜ Vh⟩ 5. A posteriori error estimation for unsteady Smagorinsky model 131 Thanks to the Lemma 5.1, and taking into account the definition of γN(µ) and βN(µ) in (5.19) and (5.20) respectively, we obtain, ⟨R(˜ UN(µ);µ),˜ Vh⟩≤ρT∥s(˜ UN(µ)−˜ Uh(µ))∥˜ X∥˜ Eh(µ)∥˜ X∥˜ Vh∥˜ X +γN(µ)∥˜ Eh(µ)∥˜ X∥˜ Vh∥˜ X. Dividing by ∥˜ Vh∥˜ X , taking supremum in ˜ Vh∈˜ Xh and using the definition of (5.31) , it holds εN(µ)≤ρT∥˜ Eh(µ)∥2 ˜ X+γN(µ)∥˜ Eh(µ)∥˜ X. Since 0 ≤τN(µ)≤1 and 1−p1−τN(µ)≤τN(µ), we have that 2ρT βN(µ)∆N(µ)≤τN(µ) and then ∆N(µ)≤2εN(µ) βN(µ). It follows that ∆N(µ)≤2ρT βN(µ)∥˜ Eh(µ)∥2 ˜ X+2γN(µ) βN(µ)∥˜ Eh(µ)∥˜ X Thanks to (5.32), we know that ∥˜ Eh(µ)∥˜ X≤∆N(µ), then 2ρT βN(µ)∥˜ Eh(µ)∥˜ X≤τN(µ). It follows that ∆N(µ)≤2γN(µ) βN(µ)+τN(µ)∥˜ Eh(µ)∥˜ X. 5.4 Stabiliy Factor approximation for unsteady problems In Theorems 5.1-5.2, we have obtained the a posteriori error bound estimator ∆N(µ) to apply the RB method through the Greedy Algorithm. Unfortunately, the construction of the a posteriori error bound estimator is expensive, specifically, due to the computation of the Stability Factor. We outline the application of the same technique as in Section 1.5 in [ 14 ] for the computation of βh(µ) for the unsteady case hereunder. 132 5.4. Stabiliy Factor approximation for unsteady problems Let us define the Stability Factor βh(µ) for the full-order problem (5.1) by βh(µ):=inf ˜ Zh∈˜ Xh sup ˜ Vh∈˜ Xh ∂1˜ A(˜ Uh(µ),˜ Vh;µ)( ˜ Zh) ∥˜ Zh∥˜ X∥˜ Vh∥˜ X (5.35) and the so called supremizer operator Tµ h:˜ Xh−→ ˜ Xh by Tµ h˜ Zh=arg sup ˜ Vh∈˜ Xh ∂1˜ A(˜ Uh(µ),˜ Vh;µ)( ˜ Zh) ∥˜ Vh∥˜ X (5.36) where ˜ Uh(µ)=(˜ uh(µ),˜ph(µ)) is the solution of problem (5.1) . Then, βh(µ) = inf ˜ Zh∈˜ Xh ∥Tµ h˜ Zh∥˜ X ∥˜ Zh∥˜ X .(5.37) We recall that the ˜ X -norm was introduced in (5.4) and it is associated to the inner product (˜ Vh,˜ Zh)˜ X= (˜ vh,˜ zh)L2 T(H1 0)+( ˜qh,˜rh)L2(QT).(5.38) where ˜ Vh= (˜ vh,˜qh) and ˜ Zh= (˜ zh,˜rh) . The key is to express (5.37) in an algebraic form and transform the computation of βh(µ) into an eigenvalue problem. This is, we look after two matrices ˜ X,˜ F(µ)∈MNhL×NhL(R) such that ˜ Vh T˜ X˜ Zh= ( ˜ Vh,˜ Zh)˜ X,∀˜ Vh,˜ Zh∈˜ Xh, ˜ Vh T˜ F(µ)˜ Zh=∂1˜ A(˜ Uh(µ),˜ Vh;µ)( ˜ Zh),∀˜ Vh,˜ Zh∈˜ Xh,∀µ∈D, where Nh=dim(˜ Xh) and vectors ˜ Vh,˜ Zh∈RNhL are the values of the degrees of freedom of the representation of ˜ Vh and ˜ Zh on ˜ Xh . For all k=1,...,L let us define Vhk= vk h qk h!,Zhk= zk h rk h!, and ˜ Vh=     V1 h V2 h . . . VL h      ,˜ Zh=     Z1 h Z2 h . . . ZL h      . Now, we rewrite (5.37) by βh(µ) = inf ˜ Zh∈˜ Xh ˜ Zh T˜ F(µ)T˜ X−1˜ F(µ)˜ Zh ˜ Zh T˜ X˜ Zh .(5.39) 5. A posteriori error estimation for unsteady Smagorinsky model 133 This quotient is the minimum of a Rayleigh quotient. Then, for a given µ∈D , we solve the following eigenvalue problem (Find (α,˜ Zh)∈Rט Xh,˜ Zh=0,such that ˜ F(µ)T˜ X−1˜ F(µ)˜ Zh=α˜ X˜ Zh,∀˜ Zh∈˜ Xh.(5.40) and from (5.39) , we obtain that βh(µ) = (αmin)1/2 with αmin the minimum eigenvalue of (5.40) . To solve the eigenvalue problem (5.40) , we perform a Power Iteration method, consequently, we shall compute the largest eigenvalue of matrix ˜ F(µ)−1˜ X˜ F(µ)−T˜ X . We know that (˜ Vh,˜ Zh)˜ X=ZTf 0(˜ vh,˜ zh)T+( ˜qh,˜rh)L2(Ω)dt = L ∑ k=1 ∆t(vk h,zk h)T+∆t(qk h,rk h)L2(Ω) Let us define the matrices L and T by (qh,rh)L2(Ω)=qhTLrh,(vh,zh)T=vhTTzh,(5.41) for all vh,zh∈Yh and qh,rh∈Qh . Then (˜ Vh,˜ Zh)˜ X= L ∑ k=1 ∆t(vk h)TTzk h+∆t(qk h)TLrk h. Then, (˜ Vh,˜ Zh)˜ X= L ∑ k=1 ∆t(Vk h)TT0 0LZk h. Finally, we obtain that (˜ Vh,˜ Zh)˜ X=˜ Vh T˜ X˜ Zh , where ˜ X=∆t       T L... T L        .(5.42) To obtain the ˜ F(µ) matrix, we recall that we have obtained the decomposition (5.15) , ∂1˜ A(˜ Uh(µ),˜ Vh;µ)( ˜ Zh) = L ∑ k=1 ∆t∂1A(Uk h(µ),Vk h;µ)(Zk h). Then, we are able to define Fk(µ) for k=1,...,L , such that for all ˜ Vh,˜ Zh∈˜ Xh and µ∈D it holds (Vk h)TFk(µ)Zk h=∂1A(Uk h(µ),Vk h;µ)(Zk h), 134 5.4. Stabiliy Factor approximation for unsteady problems and we finally define ˜ F(µ) by ˜ F(µ) = ∆t     F1(µ) F2(µ)... FL(µ)      .(5.43) We observe that matrix ˜ X in (5.42) is a block-diagonal matrix, whose elements are composed by the matrices L and T defined in (5.41) and also ˜ F(µ) whose elements are the matrices Fk(µ) for k=1,...,L . Our goal is to compute the largest eigenvalue of the matrix ˜ F(µ)−1˜ X˜ F(µ)−T˜ X which turns out to be a matrix with very large dimension. The computation of the largest eigenvalue of this matrix does not seen to be affordable with the current techniques in the literature. To sum up, the problem of building low cost techniques for the approximation of βh(µ) remains open. We intend to use hyper-reduction techniques to afford the computation of the eigenvalue problem (5.40) . We could use only some time steps in the building of ˜ F(µ) or/and we could use a coarser mesh in order to decrease Nh . 6 A posteriori error estimation based upon Kolmogórov turbulence theory The aim of this chapter is to introduce an alternative a posteriori error estimation. This estimator is based upon the Kolmogórov turbulence theory, which introduces the idea of energy cascade and an expression for the energy spectrum. In Section 6.1, we introduce an energy balance analysis of the Navier-Stokes equations to see how the energy is dissipated. In Section 6.2, we introduce the Fourier Transform and thanks to which, we are able to study the eddy energy by its size. In sections 6.3 and 6.4, we introduce the idea behind energy cascades and the energy spectrum. Finally, in Section 6.5, we introduce an a posteriori error estimator and we develop an academic test for a RB problem in Section 6.6. 6.1 System energy study In this section, we study the basics of the energy balance of the Navier-Stokes (NS) equations for three cases. Let Ω be a bounded polyhedral domain in Rd , with d=2,3 and Tf>0 be a chosen final finite time. We denote the time interval as If= (0,Tf) and the time-space domain as QT=If×Ω . The NS equations in differential form are: (∂tu+(u·∇)u=−∇p+ν∆u+f,in QT, ∇·u=0,in QT.(6.1) 135 136 6.1. System energy study We complement this problem with homogeneous or periodic boundary conditions and u0=u0(x) as the initial condition. The NS equations determine any flow by the velocity field ( [u] = ms−1 ) and pressure per mass density ( [p] = m2s−2 ). This flow can be seen as a laminar or turbulent flow, depending on the Reynolds number associated with the problem, see (A.22) . The Reynolds number establishes a relation between inertial and viscous forces, understanding that we have a turbulent flow hen inertial forces are important and laminar when viscous forces are relevant. That said, since the velocity and characteristic length are difference along the domain, we can find different Reynolds numbers within the domain. Now, we define the total kinetic energy of the system as E=1 2⟨|u|2⟩ where we denote ⟨·⟩ a suitable average operator. This average verifies the following properties: •⟨u+v⟩=⟨u⟩+⟨v⟩, •⟨au⟩=a⟨u⟩for all constant a, •⟨∂tu⟩=∂t⟨u⟩. From now on, for simplicity, we will suppose that |Ω|=1 , and we choose a spacial average as follows ⟨u⟩=ZΩ udx. We study the energy of equations (6.1) in difference situations, depending on the values of f and ν . Case 1: f= 0 and ν=0 We multiply (6.1) by u , and applying the average, we obtain that 1 2∂t⟨|u|2⟩=∂tE=0. This is, the energy is time-constant since we are removing the source and viscosity terms, indeed, E=1 2⟨|u0|2⟩ . 6. A posteriori error estimation based upon Kolmogórov turbulence theory 137 Case 2: f=0 and ν=0 We repeat the process and we obtain that ∂tE=−εν where εν=2ν⟨|Du|2⟩ is the turbulent dissipation, that we will introduce with more detail below, with Du= (1/2)(∇u+∇ut) . Note that by definition, εν≥0 , then, ∂tE≤0 . This is, the energy is decreasing over time, then we can deduce that the energy dissipation is caused by the viscosity forces. Case 3: f= 0 and ν=0 We obtain that ∂tE=εf−εν where εf=⟨f·u⟩ . In this case, we can not ensure the behavior of energy since the force terms could add or remove energy to the system, depending on the sign of εf . 6.2 Fourier space representation One of the most notorious features in turbulence flow is the variability of eddies induced by the fluid motion. Each eddy could be classified by its size r , and also by its wavenumber k=1/r . In this section, we express the energy and the turbulent dissipation over the wavenumber. This will allow us to describe the energy cascade and energy spectrum concepts in the following sections. Thanks to the Fourier transform, we are able to decompose the velocity field in Fourier modes ˆ u(t,k) as follows ˆ u(t,k) = ZΩ u(t,x)φk(x)dx and on the contrary, u(t,x) = ZRdˆ u(t,k)φk(x)dk where φk(x) = e2πik·x , φk(x) = e−2πik·x is the complex conjugate, k=kn , n a unit vector in Zd . We assume that ˆ u(t,k) = 0 for k=0 , this is, the mean flow is zero. From now on, to simplify notation, we consider u(x) = u(t,x) and ˆ u(k) = ˆ u(t,k) .