**Numerical Computations of Vortex Formation Length in Flow Past an Elliptical Cylinder**

#### **Matthew Karlson 1, Bogdan G. Nita 2,\* and Ashwin Vaidya <sup>2</sup>**


Received: 29 June 2020; Accepted: 8 September 2020; Published: 10 September 2020

**Abstract:** We examine two dimensional properties of vortex shedding past elliptical cylinders through numerical simulations. Specifically, we investigate the vortex formation length in the Reynolds number regime 10 to 100 for elliptical bodies of aspect ratio in the range 0.4 to 1.4. Our computations reveal that in the steady flow regime, the change in the vortex length follows a linear profile with respect to the Reynolds number, while in the unsteady regime, the time averaged vortex length decreases in an exponential manner with increasing Reynolds number. The transition in profile is used to identify the critical Reynolds number which marks the bifurcation of the Karman vortex from steady symmetric to the unsteady, asymmetric configuration. Additionally, relationships between the vortex length and aspect ratio are also explored. The work presented here is an example of a module that can be used in a project based learning course on computational fluid dynamics.

**Keywords:** vortex formation length; wake; vortex shedding

#### **1. Introduction**

Vortex development in a fluid's flow is a highly nonlinear phenomenon which goes through multiple bifurcations. This topic is rarely introduced in a serious manner in an undergraduate course on fluids. However, there are simple ways to talk about vortex development by combining qualitative and quantitative approaches that go beyond simply examining classical images or the use of sophisticated particle image velocimetry (PIV) techniques. We introduce one such method of talking about the physics of vortex development in this paper, which can be used in the form of a lesson plan for a lecture or to motivate computational projects in more advanced classes in fluids. Student-centered practices such as problem-based and project-based learning (PBL) are more commonly practiced in the arts. Instructional methods related to PBL promote a more inductive approach to learning whereby generalizations and abstractions follow from first understanding specific cases [1]. This approach is in contrast to the deductive strategy taken in the sciences which is a more top-down approach and a possible cause of alienation in several students. The concept of problem-based learning began more than 30 years ago in the context of medical education. PBL has been defined as the "posing of a complex problem to students to initiate the learning process" [2] and as "experiential learning organized around the investigation and resolution of messy, real-world problems" [3]. PBL can be implemented at various scales in a course with a focus from a "teacher to student-centered education with process-oriented methods of learning" [4,5]. The recent popularity of project based learning approach in physics and engineering education is based on research indicating the effectiveness of PBL in enhancing student engagement [4,6,7]. Several recent educational papers have specifically discussed the effectiveness of computational problems in fluid dynamics and the use of software such as Comsol, among others, in improving classroom engagement [8–13]. The current paper is an attempt to present a similar example of a complex problem in fluid dynamics which is apt as a unit which can be introduced as a project.

Flow past a circular cylinder is a very well studied problem in classical fluid mechanics. While the creeping flow regime is completely understood and most related problems are approachable analytically, the inertial regime contains several unanswered questions. The evolution of flow past a cylinder is a particularly interesting and well studied problem [14]. The Reynolds number (*Re* = *UL*/*ν*, where *U* and *L* represent the characteristic velocity and length respectivel, and *ν* refers to the dynamic viscosity) has been shown to capture several critical changes in flow structure. The first of these happens around *Re* = 5 [15], where the flow transitions from the creeping flow to one with a symmetric vortex profile. A second flow bifurcation from the steady symmetric profile to a asymmetric unsteady vortex occurs at around *Re* ≈ 47 [16] which lasts until about *Re* ≈ 150. Following the notation employed in the recent literature [17], we identify these critical Reynolds numbers as *Rec*<sup>1</sup> and *Rec*2, respectively. We refer the readers to the paper by Faruquee et al. [17] who provide a thorough discussion of the various historical experimental and numerical studies on this topic. The critical values referred to above are sensitive to shape of the obstacle, aspect ratio (*AR*), blockage ratio of channel diameter to channel width, roughness of the cylinder etc. but the qualitative aspects of the various transitions are still maintained.

In this paper we are concerned with highlighting the dynamics of vortex formation, specifically through an examination of the vortex length in a flow past an elliptical obstacle. Experimental measurements of the length and width of wake vortices past cylinders have been discussed in the literature [14,16,18–24]. The vortex length, usually denoted *Lw*, is defined as "the streamwise distance between the confluence point (wake stagnation point) and the rear stagnation point of the cylinder" [17]. It is also, more commonly, defined as the distance between the rear stagnation point and the "point downstream where the velocity fluctuation level has grown to a maximum" [14] (Figure 1 provides a schematic explanation of this metric). Experiments and numerics indicate the relation between *Lw* and maximum velocity fluctuation and the base pressure (at the rear of the obstacle) to be inversely proportional [23,25,26]. Coutanceau and Bouard [27] put forth the linear equation *Lw*/*d* = 0.05 *Re* for 4.4 < *Re* < 40, which, in the steady, symmetric vortex regime, related the growth of the near wake vortex as a function of *Re*.

The effect of *AR*, which is the ratio of the semi-minor (*a*) to semi-major (*b*) axis of the obstacle, is also of interest in this study. Note that *AR* = *<sup>a</sup> <sup>b</sup>* < 1 indicates an elongated body while, *AR* > 1 suggests a squat, flat object (such as disk in 3D). In two dimensions the shape itself is identical and therefore the *AR* is reflective of the orientation of the body with respect to the oncoming flow. Hence in 2D, we can describe *AR* < 1 as indicative of a high drag configuration where the entire length of the body is perpendicular to the flow, while in the case of *AR* > 1, the body assumes a pefectly streamlined position with the length of the body parallel to the flow. Recent numerical simulations [17] reveal that *Lw* has a tendency to increase with increasing *AR* for *Rec*<sup>1</sup> < *Re* < *Rec*2. However, for *AR* < 0.4, the wake disappears. For *Re* = 40, Faruquee et al. [17] provide the following equation relating the normalized vortex formation length to AR:

$$\frac{L\_w}{d} = 2.2(AR)^2 + 0.5(AR) - 0.43\tag{1}$$

where *d* is the hydraulic diameter. This formula indicates a critical minimum *AR* of 0.34 below which no symmetric vortex forms. This critical *AR* is sensitive to the *Re* and has been shown to increase with decreasing *Re* Table 3 [17]. The calculations are in agreement with the experimental studies [27,28].

**Figure 1.** This figure shows a visualization of the primary vortex region in a flow past an ellipsoidal body. Based on past experimental work [24], the vortex length (denoted by the length of the bold red line) is measured as the distance between the rear surface of the body and the "pinch-off" point where the flow velocity vanishes.

In the rest of the paper we numerically investigate the vortex formation length for a class of ellipse shaped obstacles in a flow in the range 10 < *Re* < 100 for various *AR*s. Our investigation goes past *Rec*<sup>2</sup> into the asymmetric, periodic vortex shedding regime. Section 2, discusses the numerical method used and Section 3 elaborates the results of our computations and compares them with those in the literature.

Vortex formation in flow past cylinders has been studied extensively for many practical applications. Vortices form when fluids flow past obstacles at sufficiently high speeds and are therefore of particular interest in various branches of engineering as they can pose a threat by inducing harmful vibrations in aircrafts, buildings, and other structures. In the ensuing analysis, we determine the length of the primary vortex region in the wake of an immersed body, since this region is most influential in the dynamics of the body. More specifically, we examine 2D elliptical bodies of different aspect ratio (henceforth denoted *AR*), which is a ratio of semi-major to semi-minor axis (oblate) and semi-minor to semi-major (prolate). A total of thirteen different values of *AR* from 0.2 to 2.6, increasing in increments of 0.2, have been considered. This simple numerical approach is a continuation of a previous experiment conducted on a flow past a fixed cylinder in a flow tank, in which the vortex length was determined by means of visualization and serves to elucidate a complex problem by simple means, which we believe makes for an effective class project.

#### **2. Methodology**

We used COMSOL Multiphysics to model a 2D flow in a channel past a fixed cylinder. The software uses a finite element method to solve the Navier–Stokes and incompressibility equations given by

$$
\rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) - \mu \nabla \cdot (\nabla \mathbf{u} + \nabla^T \mathbf{u}) + \nabla p = 0 \tag{2}
$$

$$
\nabla \cdot \mathbf{u} = 0.\tag{3}
$$

Here, *ρ* is the density of the fluid, **u** = (*ux*, *uy*) is the divergence free flow field, *t* is time, *μ* is the kinematic viscosity, and *p* is the pressure. The Reynolds number is defined as *Re* = *Ud <sup>ν</sup>* , where *U* is the far-field of free-stream velocity, ranging between 0.1–1 m/s, *d* is the characteristic length, which in this case was 10−<sup>1</sup> m and *ν* is the dynamic viscosity, taken to be 10−<sup>3</sup> kg/m3. The solution to the above flow equation yielded the velocity and pressure fields which were then utilized to identify the vortex length for 10 < *Re* < 100 and ellipses of varying aspect ratios (*AR* is the ratio of minor axis to major axis). The major and minor axes of the ellipses were chosen to conform to the desired values of *AR* but with keeping the area of each ellipse the same as that of the cylinder. Table 1 summarizes the important parameters used in the numerical computations.

**Table 1.** The table provides some important parameters for the numerical study. Note that the area of the elliptical cylinders is always maintained at 0.0785 m2 in all calculations.


The problem of flow past a 2D cylinder is a well studied and benchmarked problem in COMSOL [29]. A description of the code and the methodology of this standard problem can also be found on the COMSOL website: comsol.com/model/download/449401/models.mph.cylinder\_ flow.pdf. For the purposes of this work, this code was suitably adapted for the study of the elliptical cylinder. The problem was solved using the FSI module in COMSOL which uses a PARADISO solver; it was run for 5 s in increments of 0.01 s using a "fine" mesh consisting of 5662 elements. Convergence tests were performed in a previous study [30,31] for a more complex problem involving 2D and 3D elliptical cylinders with attached wings (or flexible fiber).

The dimensions and geometry of the problem studied here were similar to those used in our previous studies [30]. We also verified our code for the case of flow past a circular cylinder with perfect slip conditions along the top and bottom walls and no slip on the surface of the cylinder, to mimic unbounded flow. Through this calculation we were able to obtain the critical bifurcation of around 47 where the flow transitions from steady to unsteady. This helped us to confirm that the metric used to identify the critical Reynolds number was correct. In a classroom setting the relationship between the Reynolds number and the Strouhal number would be appropriate as an additional validation of the numerical scheme. The remaining computations were conducted in a bounded channel with no-slip conditions on the channel walls and also on the surface of the ellipse. The various panels in Figures 2 and 3 show some sample flow and pressure profiles based on our numerical simulations for various Reynolds numbers.

Two different criteria were used to measure the length of the vortex depending on whether the flow is steady or unsteady. For the steady case, since the flow does not vary in time (see the top panels in Figure 4), the vortex region is symmetric about *y* = 0 and its length remains constant in time. In this case, we examined the horizontal component of the velocity, *ux*, and of the pressure gradient, namely, *dp dx* as a function of *x* along *y* = 0. The border or separatrix of the primary vortex pair separates the inner wake flow from the outer uniform flow. At the junction of the separatrix and *y* = 0, the two flows are in opposing directions, the outer flow moving in the positive *x* while the wake flow moves in the negative *x* direction. The point at which *ux* = 0 can therefore be defined as a critical point which marks the end of the primary vortex region. The horizontal distance of this critical point from the rear end of the cylinder is defined as the vortex length. To confirm this hypothesis, a second criterion was also tested. We also observe that the pressure along *y* = 0 increases the fastest as we cross the critical point. Therefore the vortex length can also be defined with respect to the maximum of the magnitude of the pressure gradient. Since *dp dy* is negligible, this collapses to the maximum of *dp dx* along *y* = 0.

(**a**) Streamlines of the flow in the steady regime (**b**) Pressure profile in the wake of the ellipse

**Figure 2.** Numerical simulations of the flow past elliptical objects of various *AR*s: steady flow.

(**a**) Streamlines of the flow in the unsteady regime (**b**) Pressure profile in the wake of the ellipse

**Figure 3.** Numerical simulations of the flow past elliptical objects of various *AR*s: unsteady flow.

**Figure 4.** The panels (**a**,**b**) show the position of the wake stagnation point in the case of steady flow, while panels (**c**,**d**) show the periodically changing position of the stagnation point for the oscillating flow regime.

In the case when flow becomes unsteady, the wake region oscillates in space and time (see the bottom panels in Figure 4). The vortex region is asymmetric about *y* = 0 and moves in both spatial directions *x* and *y*. In this case the critical point at a given time was defined by the maximum of the magnitude of the pressure gradient, ∇*p*, which is a generalization of the criterion used in the steady state case and the appropriate length is the time averaged vortex length over a period. We use numerical interpolation to smooth the data and find our zeros and extremum.

#### **3. Results and Discussion**

The computations for the steady case show a monotonic increase in vortex length with increasing *Re*, below a critical threshold, for all *AR*s (see for example, Figure 5). The rising trend is similar for both criteria employed to measure the vortex length. When the flow becomes unsteady, the mean vortex length is seen to decrease with increasing *Re*. The maximum vortex length can be associated with the onset of a transition in stability as the flow bifurcates from steady to unsteady Karman vortex. The Figure 5a–d show the mean change in vortex length versus *Re* for some sample cases of flow past ellipses of *AR* 0.4, 0.8, 1.0, and 1.2. All graphs, show the same overall profile and the turning point can be used to identify the critical Reynolds number, *Rec* which is sensitive to the *AR* of the body. The Table 2 depicts the *Lw*/*d* versus *Re* profile for 10 < *Re* < 60, i.e., for flows that generate a steady vortex, which is in keeping with the linear relationship proposed by Coutanceau and Bouard [27]. The table shows the slopes of best fit lines for various *AR*s which appear to monotonically decrease with increasing *AR* and are overall close to the value suggested in the literature. In the regime 60 < *Re* < 100, we use an exponential decaying function, *Lw*/*d* = *AeαRe* to fit the data. Table 2 also shows the value of the exponent, *α*, along with the goodness of fit for various AR. As in the case of the steady vortex, the dependence of the vortex formation length on *Re* is sensitive to AR. The missing data for *AR* = 1.4 is due to the termination of our study at *Re* = 100; at this point the *AR* = 1.4 case does not reveal an unsteady vortex. Table 3 also shows the average *Lw*/*d* for different *AR* and *Re*.

**Figure 5.** The figures show the change in normalized vortex length, *Lw*/*d*, versus *Re* for *AR* = 0.4 in (**a**), *AR* = 0.8 in (**b**), *AR* = 1.0 in (**c**) and *AR* = 1.2 in (**d**). The maximum point in each graph characterizes the transition *Re* at which the flow changes from steady symmetric to unsteady asymmetric.

**Table 2.** This table shows the nature of the *Lw* vs. *Re* relationship in the two different flow regimes. The first two rows provide the slope of the best line representing the changing value of *Lw* in the steady wake. Rows three and four tabulate the exponentially decaying rate of *Lw* as a function of *Re* in the unsteady, asymmetric wake.


A particularly interesting relationship to investigate in this problem is that between the critical Reynolds numbers, *Rec* and *AR*. The transitions from steady to unsteady wake are very sensitive to the specific geometric characteristics of the cylinder. Figure 6 shows the critical *Re* for *AR*s ranging between 0.4 and 1.4 estimated from using the velocity and pressure gradient criteria. The *Rec* vs. *AR* curve shows a minimum at *AR* ≈ 0.6–0.7. Overall, for bodies with *AR* in the range 0.4–1.4, *Rec* appears to lie in the range 55–95; for a cylinder our calculations show this critical value to lie at the accepted value of *Rec* = 70. The more streamlined bodies show greater propensity to develop elongated primary vortices. The equation:

$$\text{Re}\_{crit} = 79.46(AR)^2 - 116.18(AR) + 100.26 \tag{4}$$

captures the quadratic nature of the profile (with *R*<sup>2</sup> value of 0.96) seen in Figure 6, obtained by fitting the average of the two curves shown in the figure. This equation gives us valuable information about the optimal geometry and its relationship with vortex formation. Engineers and designers wanting to suppress disruptive osculations may find this equation particularly useful.

**Figure 6.** The graph shows the variation of *Recrit* versus *AR*.

**Table 3.** The table shows the average values of the normalized vortex length *Lw*/*d* for various values of *AR* and *Re*.


#### **4. Conclusions**

In summary, our overall computations and qualitative profiles are in agreement with previous experimental results [24], where measurements of variations in vortex length were made through flow visualization and imaging techniques of flow past various cylinders in a flow tank. These experiments indicate a similar extremum in vortex length as a function of flow speed (*Re*) and *AR*. The central contribution of this work lies in two consistent ways of defining the vortex formation length and the observation of the vortex length dependence on *Re*, especially in the unsteady flow regime. Questions for future investigation include extending the *Re* regime of the study and experimental verification of the unsteady vortex formation length. Pedagogical treatments in the literature appear to focus primarily on lower level courses in science and engineering. There is little discussion about the teaching and learning of more advanced topics like fluid dynamics [32]. We believe that examples like the ones presented here add value to the science of fluid dynamics in some capacity, but also to the ways in which students can engage with difficult topics in fluid dynamics. Such a treatment lends itself very well to a first serious course in fluids but also in more advanced courses where students have greater proficiency in dealing with computational software like Comsol.

**Author Contributions:** The computations in this paper were performed while M.K. was a graduate student at Montclair State University. Authors B.G.N. and A.V. designed the study, helped analyze the results and wrote the paper. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## *Article* **On the Existence of Leray-Hopf Weak Solutions to the Navier-Stokes Equations**

**Luigi C. Berselli 1,\* and Stefano Spirito <sup>2</sup>**

	- <sup>2</sup> DISIM—Dipartimento di Ingegneria e Scienze dell'Informazione e Matematica,

**Abstract:** We give a rather short and self-contained presentation of the global existence for Leray-Hopf weak solutions to the three dimensional incompressible Navier-Stokes equations, with constant density. We give a unified treatment in terms of the domains and the relative boundary conditions and in terms of the approximation methods. More precisely, we consider the case of the whole space, the flat torus, and the case of a general bounded domain with a smooth boundary (the latter supplemented with homogeneous Dirichlet conditions). We consider as approximation schemes the Leray approximation method, the Faedo-Galerkin method, the semi-discretization in time and the approximation by adding a Smagorinsky-Ladyžhenskaya term. We mainly focus on developing a unified treatment especially in the compactness argument needed to show that approximations converge to the weak solutions.

**Keywords:** Navier-Stokes equations; Leray-Hopf weak solutions; existence

**Citation:** Berselli, L.C.; Spirito, S. On the Existence of Leray-Hopf Weak Solutions to the Navier-Stokes Equations. *Fluids* **2021**, *6*, 42. https://doi.org/10.3390/fluids 6010042

Received: 15 December 2020 Accepted: 5 January 2021 Published: 13 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Let *<sup>T</sup>* <sup>&</sup>gt; 0 be an arbitrary finite number representing the time, <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup><sup>3</sup> be a domain to be specified later, and *ν* > 0 be a positive number representing the kinematic viscosity. The incompressible Navier-Stokes equations model the dynamic of a viscous and incompressible fluid at constant temperature and with constant density. They are given by the following system of PDE's posed in (0, *T*) × Ω:

$$\begin{cases} \partial\_t u + \left( u \cdot \nabla \right) u + \nabla p - \nu \Delta u = f & \text{in } (0, T) \times \Omega, \\ \qquad \text{div} \, u = 0 & \text{in } (0, T) \times \Omega. \end{cases} \tag{1}$$

The vector field *<sup>u</sup>* <sup>∈</sup> <sup>R</sup><sup>3</sup> is the velocity, *<sup>p</sup>* <sup>∈</sup> <sup>R</sup> is the scalar pressure, and to avoid inessential complications, we set the external force *f* = 0 (but all results presented here can be easily extended to the case of a non vanishing external force, see Remark 4). The first equation is the conservation of linear momentum and the second equation, also called the incompressibility constraint, can be considered as the conservation of the mass, since the density is assumed to be constant. The system (1) has to be supplemented with initial and boundary conditions. Regarding the initial condition we impose that

$$u|\_{t=0} = u\_{0\prime} \quad \text{in } \Omega\_{\prime}$$

with *u*<sup>0</sup> satisfying the compatibility condition div *u*<sup>0</sup> = 0 in Ω. For the boundary conditions we need to specify the assumptions on the domain. We consider three cases, Ω = R3, <sup>Ω</sup> <sup>=</sup> <sup>T</sup><sup>3</sup> with <sup>T</sup><sup>3</sup> being the three-dimensional flat torus, and <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup><sup>3</sup> being a bounded domain, whose boundary will be denoted by *∂*Ω; we refer to Assumption 1 for the precise hypotheses on Ω.

For each of the three different cases we impose the different and natural boundary conditions:

(*i*) *<sup>u</sup>* <sup>→</sup> 0 for <sup>|</sup>*x*| → <sup>∞</sup>, if <sup>Ω</sup> <sup>=</sup> <sup>R</sup>3;

$$\mu(ii) \quad \mu \text{ is periodic, if } \Omega = \mathbb{T}^3; \tag{2}$$

(*iii*) *u* = 0 on (0, *T*) × *∂*Ω, if Ω is a bounded domain.

Note that the initial datum will be requested to be tangential to the boundary in the case (*i*), and to satisfy the condition (*ii*) and (*iii*) in the other cases. Contrary to the system of compressible Navier-Stokes equations, the pressure *p*, instead of being obtained through a state equation, is an unknown of the system. This is a consequence of the incompressibility conditions and indeed the pressure can be interpreted as Lagrange multiplier associated with the incompressibility constraint. Note that there are no initial/boundary conditions imposed on the pressure, which (since it appears only as a gradient in the momentum equation) is always determined up to an arbitrary function of time.

Generally speaking, it is very difficult to prove existence and uniqueness of smooth solutions to nonlinear PDE's. Here, with existence we always mean global in time existence, namely existence on any given time interval (0, *T*), for arbitrary *T* > 0. The available theories for weak solutions provide a framework to give a proper meaning to PDE's, without requiring too much regularity on the solutions and they rely on the theory of generalized functions and distributions. In particular, the landmark idea in the theory of weak solutions is to give up on solving the equations point-wise but trying to solve them in an averaged sense, which is meaningful also from a physical point of view. In the case of fluid mechanics, we expect a very complex behavior by (turbulent) flows appearing in real life, hence we expect to be able to capture only averages of the velocity and pressure, see Reference [1].

The problem of global existence generally becomes easier since the class of available solution is enlarged and several functional analysis tools can be now used. However, the price to pay to have such a relatively simple existence theory is that the uniqueness problem becomes a very difficult one and many calculation which are obvious when dealing with smooth solutions are not possible or hard to be justified. The three-dimensional incompressible Navier-Stokes is a paradigmatic example of a such situation and the introduction of weak solutions dates back about 100 years ago. In fact, in a series of celebrated papers (the time evolution is treated in Reference [2]) Jean Leray introduced the notion of weak solution as a mathematical tool, but also with a strong understanding of the physics behind the equations. The theory of weak solutions is also strictly linked with the name of Eberhard Hopf [3] who gave the first contribution to the problem of existence of weak solutions in a bounded domain, by means of the Faedo-Galerkin method.

It is interesting to observe that many methods and techniques of functional analysis (which are now a common background of graduate students in mathematics) originated from the study of PDE's and especially from those arising in fluid mechanics. In this note we are trying to explain an extremely limited part of the theory: the existence of (Leray-Hopf) weak solutions. This is a topic at the level of most undergraduate students, with a minimal knowledge of Sobolev spaces and functional analysis (mainly weak convergence and weak compactness), as for instance in the widely used (text)books by Brezis [4], just to name one. Note also that we try to present a minimal spot in the abstract theory of Navier-Stokes equations, which can be an "appetizer" for students trying to start a serious understanding of (part of) the mathematical fluid mechanics. It is impossible to review what is done on the subject, even only for the mathematical analysis side. Nevertheless many information, at an introductory or more advanced level, can be found in several books, see for instance, just to name a few in alphabetical order [5–11].

We think that we will not discourage any reader unfolding the (many) mathematical difficulties of the topic, but –instead– we hope that highlighting the challenges which are typical of mathematical fluid mechanics further interest could be stimulated; To this end we quote the following coming from an interview reported in Reference [12] in a essay in memory of Jacob Schwartz:

*When I asked him [Jacob Schwartz] if there was a subject he had trouble learning, he admitted that there was, namely, fluid dynamics. "It is not a subject that can be expressed in terms of theorems and their proofs," he said.*

Following the above point of view, the first step even in the mathematical analysis of the Navier-Stokes equations is that of giving an appropriate definition of weak solutions, which take into account the functional spaces where it is reasonable find weak solutions, the initial and the boundary conditions. Usually, the functional space to be considered are hinted by the a priori estimate available for the system under consideration. The informal notion of an a priori estimate may be a quantitative bound depending only on the data of the problem, which holds for smooth solutions of the system under consideration, regardless their existence. In particular, for system arising from physics, the a priori estimates usually have a deep physical interpretation.

In the context of the three-dimensional incompressible Navier-Stokes equations the main a priori estimate is indeed the conservation of the energy of the system and is given by the following integral equality:

$$\int |u(t, \mathbf{x})|^2 \, d\mathbf{x} + 2\nu \int\_0^t \int |\nabla u(s, \mathbf{x})|^2 \, d\mathbf{x} ds = \int |u\_0(\mathbf{x})|^2 \, d\mathbf{x} \qquad t \in [0, T], \tag{3}$$

where the space integral is over the domain under consideration.

The equality (3) has a very simple *formal* proof. Indeed, let (*u*, *p*) be a smooth solution of (1) and (2). By multiplying the momentum equation by *u* and integrating over Ω we get

$$\int \partial\_t \boldsymbol{u} \cdot \boldsymbol{u} - \nu \Delta \boldsymbol{u} \cdot \boldsymbol{u} + \left(\boldsymbol{u} \cdot \nabla\right) \boldsymbol{u} \cdot \boldsymbol{u} + \nabla p \cdot \boldsymbol{u} \, d\boldsymbol{x} = 0.$$

By integrating by parts and using the divergence free condition and (2) we get

$$
\partial\_t - \int \Delta u \cdot \boldsymbol{u} \, d\mathbf{x} = \int |\nabla \cdot \boldsymbol{u}|^2 \, d\mathbf{x}, \qquad \int (\boldsymbol{u} \cdot \nabla) \, \boldsymbol{u} \cdot \boldsymbol{u} \, d\mathbf{x} = 0, \qquad \text{and} \qquad \int \nabla p \cdot \boldsymbol{u} \, d\mathbf{x} = 0.
$$

Then, after integration in time on (0, *t*) with *t* ∈ (0, *T*) we get (3). Note that (3) gives a quantitative bound depending only on *T*, and *u*<sup>0</sup> of square integrals of the velocity field *u* and its gradient ∇*u*. The energy equality (3) will serve as motivation for the definition of Leray-Hopf weak solution we will give in Section 3.

Once a reasonable definition of weak solution is given, to prove global existence one usually exploits what it is know as a *compactness argument*, which consists in (1) proving the existence of a sequence of relatively smooth approximating solutions satisfying appropriate uniform estimates; (2) proving that limits of these approximating solutions are effectively weak solution of the problem under consideration. We remark that usually the uniform bounds obtained on the sequence of approximating solutions are the same inferred by the a priori estimates available for the system under consideration; These bounds are then hopefully inherited by weak solutions obtained with a passage to the limit. To be more precise, in the case of the Navier-Stokes equations, the approximation method should be chosen such that the approximate solutions satisfy the energy (in)equality. Due to the limited regularity which can be generally inferred on weak solutions, the validity of any energy balance on the weak solutions to the 3D Navier-Stokes equations is obtained with a limiting process on the approximate solutions and not using the solution *u* itself as a test function as done to obtain (3), since this argument is only formal and not justified when dealing with genuine Leray-Hopf weak solutions.

In this short note we provide a rather self-contained account on the global existence of weak solutions for the three-dimensional incompressible Navier-Stokes equations and some of the (several) approximation methods used in the literature. Since the convergence argument is essentially the same for every approximation methods and for every choice of the domains and boundary conditions mentioned above, we introduce (for the purpose of the exposition) a notion of *approximating solution* for which we will prove the convergence

to a Leray-Hopf weak solution of the problem (1) and (2). This is not the historical path, but is a way we identify to have a unified treatment, which can describe the existence theory within the notion of *approximating solutions.*

Then, we show how several and well-known approximations fit in the framework introduced and, therefore, we recover the existence of Leray-Hopf weak solution by using those methods. In particular, we will consider the most common techniques available for the construction. Further results based on the energy type methods, concerning uniqueness, regularity and the connection with applied analysis of turbulent flows, can be found in the forthcoming monograph [1], which is also written in the spirit of being an introduction for undergraduate students, interested in applied analysis of the Navier-Stokes equations.

#### *Organization of the Paper*

The paper is organized as follows: In Section 2 we introduce the functional spaces that we use. Then, in Section 3 we define of Leray-Hopf weak solutions and study their main properties. In Section 4 we give the definition of *approximating solution* and we prove the convergence to a Leray-Hopf weak solution. Finally, in Section 5 we prove that certain approximating schemes fit in the framework of *approximating solution*.

#### **2. Preliminaries**

In this section we fix some notations and we recall some basic preliminaries we will need for the analysis. We start by fixing the assumptions on the domain Ω.

**Assumption 1.** *The domain* <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup><sup>3</sup> *will be of the following type:*


#### *2.1. Notation*

We will never distinguish between scalar and vector functions unless it is not clear from the context. We will denote by *C*<sup>∞</sup> *<sup>c</sup>* (Ω) the space of compactly supported functions which are infinitely differentiable and D- (Ω) its dual, which is the space of distributions over Ω. In the case Ω = T<sup>3</sup> the subscript "*c*" is not needed and we set *C*<sup>∞</sup> *<sup>c</sup>* (T3) = *<sup>C</sup>*∞(T3). With an abuse of notation we will use *C*<sup>∞</sup> *<sup>c</sup>* (Ω) for all the three choices of the domain Ω satisfying Assumption 1. We recall that for any vector *<sup>f</sup>* <sup>∈</sup> *<sup>C</sup>*<sup>∞</sup> *<sup>c</sup>* (Ω; R3) the Helmholtz decomposition holds true: there exists two function *<sup>g</sup>* <sup>∈</sup> *<sup>C</sup>*<sup>∞</sup> *<sup>c</sup>* (Ω; <sup>R</sup>3) and *<sup>q</sup>* <sup>∈</sup> *<sup>C</sup>*<sup>∞</sup> *<sup>c</sup>* (Ω; R) such that *f* = *g* + ∇*q*, and *g* is divergence-free. Given a Banach space *E*, we denote with · *<sup>E</sup>* its norm. However, for the classical Lebesgue spaces *<sup>L</sup>p*(Ω), with *<sup>p</sup>* <sup>∈</sup> [1, <sup>∞</sup>], we shall denote their norms with · *p*. Finally, we recall that the space *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω) is the classical Sobolev space obtained as a closure of *C*<sup>∞</sup> *<sup>c</sup>* (Ω) in the norm

$$\|v\|\_{H^1} := \left(\int\_{\Omega} |v|^2 + |\nabla v|^2 \, dx\right)^{\frac{1}{2}}, \qquad v: \Omega \mapsto \mathbb{R}^k.$$

The subscript "0" is needed only when Ω is a bounded domain. In the case of Ω = R<sup>3</sup> or Ω = T<sup>3</sup> we have *H*<sup>1</sup> <sup>0</sup> (R3) = *<sup>H</sup>*1(R3) and *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (T3) = *<sup>H</sup>*1(T3), but as before, with an abuse of notation, we will use *H*<sup>1</sup> <sup>0</sup> (Ω) for each one of the three choices of the domain Ω satisfying Assumption 1. Moreover, we recall that *H*1(R3) and *H*1(T3) can also be characterized in terms of the Fourier Transform and the Fourier Series, respectively. When dealing with a Banach space (*E*, . *<sup>E</sup>*) we denote by *xn* → *x*, *xn x* and *xn* <sup>∗</sup> *x*, the strong, weak and weak\* convergence, respectively.

Next, let *<sup>E</sup>* be a Banach space, then *<sup>L</sup>p*(0, *<sup>T</sup>*; *<sup>E</sup>*), with 1 <sup>≤</sup> *<sup>p</sup>* <sup>&</sup>lt; <sup>∞</sup>, and *<sup>L</sup>*∞(0, *<sup>T</sup>*; *<sup>E</sup>*) denote the classical Bochner spaces of strongly measurable (classes of) functions *u* : (0, *T*) → *E* such that

$$\begin{aligned} \|\mu\|\_{L^p(E)} &:= \left(\int\_0^T \|\mu(s)\|\_E^p \, ds\right)^{\frac{1}{p}} < \infty, \\\|\mu\|\_{L^\infty(E)} &:= \operatorname\*{ess\,sup}\_{t \in [0,T]} \|\mu(t)\|\_E < \infty. \end{aligned}$$

Finally, the space of weakly continuous functions in *E*, which is denoted by *Cw*([0, *T*]; *E*), consists of functions *u* : [0, *T*] → *E* such that for any *f* ∈ *E*<sup>∗</sup> the real function of real variable

$$\langle f, \mathfrak{u} \rangle\_{E^\* \times E} : \ [0, T] \ni \mathfrak{t} \longmapsto \langle f, \mathfrak{u}(\mathfrak{t}) \rangle\_{E^\* \times E \times E}$$

is continuous.

Finally, when we write *A B*, this means that there exists a constant *c* > 0 (independent on the relevant parameters of the problem) such that *A* ≤ *cB*.

#### *2.2. The Spaces H and V*

In the analysis of solutions of the Navier-Stokes equations is useful to consider spaces of divergence-free functions. We start by defining the space

$$\mathcal{V}(\Omega) := \{ \phi \in \mathcal{C}\_{\mathfrak{c}}^{\infty}(\Omega) : \text{div}\,\phi = 0 \},$$

Then, we define the spaces

$$H := \overline{\mathcal{V}(\Omega)}^{\|\cdot\|\_2}, \qquad V := \overline{\mathcal{V}(\Omega)}^{\|\cdot\|\_{1,2}}.$$

We start by noticing that *H* and *V* are closed subspace of *L*2(Ω) and *H*<sup>1</sup> <sup>0</sup> (Ω), respectively. Therefore, they are Hilbert space themselves with the inherited scalar products, which are

$$(u, v) := \int\_{\Omega} u \cdot v \, d\mathbf{x} \qquad ((u, v)) := \int\_{\Omega} u \cdot v + \nabla u : \nabla v \, d\mathbf{x}.$$

Next, although *H* and *V* are Hilbert space, hence reflexive, we will not identify them with their duals. We will instead denote by *H* and *V* the topological dual of *H* and *V* endowed with the classical dual norms

$$\|f\|\_{H'} := \sup\_{\phi \in \mathcal{V}(\Omega), \|\phi\|\_{2} \le 1} |\langle f, \phi \rangle\_{H' \times H}|\_{\ast}$$

$$\|f\|\_{V'} := \sup\_{\phi \in \mathcal{V}(\Omega), \|\phi\|\_{1,2} \le 1} |\langle f, \phi \rangle\_{V' \times V}|\_{\ast}$$

We stress that *H* and *V* are not subset of the space of distributions D- (Ω) since D(Ω) ⊂ *H*.

Finally, we recall that by Sobolev embedding theorem and the interpolation inequality for the *<sup>L</sup>p*-norm, there exists a constant *<sup>C</sup>* <sup>&</sup>gt; 0 such that for any 2 <sup>≤</sup> *<sup>p</sup>* <sup>≤</sup> 6, *<sup>θ</sup>* <sup>=</sup> <sup>6</sup>−*<sup>p</sup>* <sup>2</sup>*<sup>p</sup>* ∈ [0, 1], and any *<sup>u</sup>* <sup>∈</sup> *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω) it holds that

$$\|\|u\|\|\_{p} \le \mathbb{C} \|\|u\|\|\_{2}^{\theta} \|\nabla u\|\|\_{2}^{1-\theta}. \tag{4}$$

The inequality (4) is a particular case of the well-known Gagliardo-Nirenberg-Sobolev inequality, see Reference [13].

#### **3. Definition of Leray-Hopf Weak Solutions**

In this section we give the definition of Leray-Hopf weak solutions and we prove some related properties. The definition is the following.

**Definition 1.** *A measurable vector field <sup>u</sup>* : (0, *<sup>T</sup>*) <sup>×</sup> <sup>Ω</sup> → <sup>R</sup><sup>3</sup> *is a Leray-Hopf weak solution of the Navier-Stokes Equations* (1) *and* (2) *if the following conditions are satisfied.*

*1. It holds that*

$$\forall u \in \mathbb{C}\_{w}([0, T]; H) \cap L^{2}(0, T; V); \tag{5}$$

*2. For any <sup>φ</sup>* ∈ V(Ω) *and any <sup>χ</sup>* <sup>∈</sup> *<sup>C</sup>*∞([0, *<sup>T</sup>*))*, it holds that*

$$\int\_{0}^{T} \left( (u(t), \phi) \dot{\chi}(t) - ((u \cdot \nabla)u, \phi) \chi(t) - \nu(\nabla u, \nabla \phi) \chi(t) \right) dt - (u\_0, \phi) \chi(0) = 0; \tag{6}$$

*3. For any t* ∈ [0, *T*]

$$\|\|u(t)\|\|\_{2}^{2} + 2\nu \int\_{0}^{t} \|\nabla u(s)\|\_{2}^{2} \, ds \le \|\|u\_{0}\|\|\_{2}^{2}. \tag{7}$$

**Remark 1.** *It is important to point out that it is an open problem whether or not condition* (7) *can be deduced from the conditions* (5) *and* (6)*. Note also that in the definition we have* (7) *which is the so-called global energy inequality and not the equality* (3)*.*

**Remark 2.** *In literature Leray-Hopf weak solutions are often defined in the space L*∞(0, *T*; *H*) *rather than Cw*([0, *T*]; *H*) *and satisfying* (7) *for a.e. everywhere t* ∈ (0, *T*) *instead that for any t* ∈ (0, *T*)*. This is equivalent to Definition 1, because in that case the velocity field can redefined on a set of measure zero in time in order to lie in Cw*([0, *T*]; *H*) *and satisfying* (7) *for any t* ∈ (0, *T*)*, see Reference [14]. We preferred to start with a solution already weakly continuous, to avoid the technical step of redefinition.*

We want to show that once we have proved the existence of a vector field satisfying the conditions in the Definition 1, we are actually solving the initial value boundary problem (1) and (2) in the sense of distributions. First of all we notice that from the condition (1), we can deduce that *u* is divergence-free and satisfies the boundary conditions (2) in the appropriate weak sense. The following lemma guarantee that *u* attains the initial datum *u*0.

**Lemma 1.** *Let u*<sup>0</sup> ∈ *H and u a Leray-Hopf weak solution. Then,*

*u*(*t*) → *u*<sup>0</sup> *strongly in H*.

**Proof.** For *<sup>k</sup>* <sup>∈</sup> <sup>N</sup> and ¯*<sup>t</sup>* <sup>∈</sup> (0, *<sup>T</sup>*), we consider the following function

$$\chi\_k^{\vec{t}}(t) = \begin{cases} \begin{array}{c} 1, & t \in [0, \vec{t}) \\ k(\vec{t} - t) + 1, & t \in [\vec{t}, \vec{t} + \frac{1}{k}) \\ 0, & t \in [\vec{t} + \frac{1}{k'}T). \end{array} \end{cases}$$

Then, by using *χ*¯*<sup>t</sup> <sup>k</sup>*, after sending *k* → ∞ and using that *u* ∈ *Cw*([0, *T*]; *H*) we arrive to the following estimate:

$$\begin{split} \left| \left( u(t), \phi \right) - \left( u\_0, \phi \right) \right| &\leq \int\_0^t \left| \left( \nabla \ u(s), \phi \right) \right| + \left| \left( u(s) \cdot \nabla \right) u(s), \phi \right| \, ds \\ &\leq t^{\frac{1}{2}} \left( \int\_0^T \left\| \nabla u(s) \right\|\_{2}^2 \right)^{\frac{1}{2}} \left( \left\| \phi \right\|\_{2}^{\frac{1}{2}} + \left\| \phi \right\|\_{\infty} \sup\_{s \in \left[ 0, T \right]} \left\| u(s) \right\|\_{2} \right). \end{split}$$

Then, for any fixed *<sup>φ</sup>* ∈ V(Ω) we can send *<sup>t</sup>* <sup>→</sup> <sup>0</sup><sup>+</sup> and we can conclude that (*u*(0), *φ*)=(*u*0, *φ*). By using the Helmholtz decomposition we deduce that this is true for any *<sup>φ</sup>* <sup>∈</sup> *<sup>C</sup>*<sup>∞</sup> *<sup>c</sup>* (Ω)) and therefore *u*(0) = *u*<sup>0</sup> a.e. on Ω. Moreover, the previous calculations also show that

$$
\mu(t) \to \mu\_0 \qquad \text{in } \mathbb{C}\_w([0, T]; H) \qquad \text{as } t \to 0^+.
$$

and, again by using the Helmholtz decomposition, the same result is valid also for *φ* ∈ *L*2(Ω). By weak lower semi-continuity of norms in weak convergence we get

$$\|\mu\_0\|\_2 \le \liminf\_{t \to 0^+} \|\mu(t)\|\_2^2.$$

Next, by using the energy inequality (7) we also get, by disregarding the non-negative dissipative term and taking the superior limit that

$$\limsup\_{t \to 0^+} \|u(t)\|\_2^2 \le \|u\_0\|\_2^2.$$

This shows that *u*(*t*) <sup>2</sup> → *u*<sup>0</sup> 2, which combined with the weak convergence implies the strong convergence, since we are in an Hilbert space. Since the norm induced on *H* is the same as in *L*2(Ω), this proves the strong convergence also in *H*.

Finally, we show that to any Leray-Hopf weak solution *u* it is possible to associate a pressure *p* such that (*u*, *p*) solves the momentum equation in (1) in the sense of distributions.

**Lemma 2.** *Let u be a Leray-Hopf weak solution of* (1) *and* (2)*. Then, there exists p* ∈ D- ((0, *T*); ×Ω)) *such that*

$$
\partial\_t \left( \mu - \nu \Delta \,\mu + (\mu \cdot \nabla) \,\mu + \nabla \,p = 0 \quad \text{in } \mathcal{D}'((0,T) \times \Omega).
$$

*and, for any t* <sup>∈</sup> (0, *<sup>T</sup>*)*, we have p*(*t*) <sup>∈</sup> *<sup>L</sup>*<sup>2</sup> *loc*(Ω) *and* ( <sup>Ω</sup> *p*(*t*) *dx* = 0*.*

In the case of a general bounded domain Ω satisfying the Assumption 3, the proof of the Lemma 2 is very technical and requires several preliminaries of operator theory. We refer to References [7,10,11,15,16] for the proof. On the other hand in the case of Ω has no physical boundary the proof is straightforward. We consider here the case Ω = T3.

**Proof.** Let *u* be a Leray-Hopf weak solution in the sense of Definition 1. For a.e. *t* ∈ (0, *T*) consider the elliptic problem

$$\begin{aligned} -\Delta p(t) &= \text{div}(u(t) \cdot \nabla u(t)) \qquad \text{in } \mathbb{T}^3, \\ \int\_{\mathbb{T}^3} p(t) \, dx &= 0. \end{aligned} \tag{8}$$

Note that by (5), Gagliardo-Nirenberg Sobolev inequality (4), and standard elliptic regularity we can infer that there exists a unique solution of (8) satisfying *<sup>p</sup>* <sup>∈</sup> *<sup>L</sup>*<sup>5</sup> <sup>3</sup> ((0, *<sup>T</sup>*) <sup>×</sup> <sup>T</sup>3). Next, we show that (*u*, *p*) solve the Navier-Stokes equations in the sense of distributions. Let *<sup>ψ</sup>*(*t*, *<sup>x</sup>*) = *<sup>χ</sup>*(*t*)*φ*(*x*) with *<sup>χ</sup>* <sup>∈</sup> *<sup>C</sup>*<sup>∞</sup> *<sup>c</sup>* (0, *<sup>T</sup>*) and *<sup>φ</sup>* <sup>∈</sup> *<sup>C</sup>*∞(T3). Let *<sup>φ</sup>* <sup>=</sup> *<sup>P</sup><sup>φ</sup>* <sup>+</sup> *<sup>Q</sup><sup>φ</sup>* be the Helmholtz decomposition, where we denote by *Pφ* the divergence-free part of *φ*. Then, since *P* and *Q* commute with derivatives because there are no physical boundaries, we have that

$$\begin{aligned} &\int\_{0}^{T} \left( u(t), \boldsymbol{\phi} \right) \dot{\boldsymbol{\chi}}(t) - \left( \left( \boldsymbol{u} \cdot \nabla \right) \boldsymbol{u}, \boldsymbol{\phi} \right) \boldsymbol{\chi}(t) - \nu \left( \nabla u, \nabla \boldsymbol{\phi} \right) \boldsymbol{\chi}(t) + \left( p(t), \text{div } Q\boldsymbol{\phi} \right) \boldsymbol{\chi}(t) \, dt \\ &= \int\_{0}^{T} \left( u(t), P\boldsymbol{\phi} \right) \dot{\boldsymbol{\chi}}(t) - \left( \left( \boldsymbol{u} \cdot \nabla \right) \boldsymbol{u}, P\boldsymbol{\phi} \right) \boldsymbol{\chi}(t) - \nu \left( \nabla u, \nabla P\boldsymbol{\phi} \right) \boldsymbol{\chi}(t) \\ &- \int\_{0}^{T} \left( \left( \boldsymbol{u} \cdot \nabla \right) \boldsymbol{u}, Q\boldsymbol{\phi} \right) \boldsymbol{\chi}(t) - \left( p(t), \text{div } \boldsymbol{\phi} \right) \boldsymbol{\chi}(t) \, dt \\ &= - \int\_{0}^{T} \left( \left( \boldsymbol{u} \cdot \nabla \right) \boldsymbol{u}, Q\boldsymbol{\phi} \right) \boldsymbol{\chi}(t) - \left( p(t), \text{div } Q\boldsymbol{\phi} \right) \boldsymbol{\chi}(t) \, dt = 0, \end{aligned} \tag{9}$$

where we have used (6) in the second equality and (8) together with the fact that *Qφ* = ∇*q* for some *<sup>q</sup>* <sup>∈</sup> *<sup>C</sup>*∞(T3) in the last equality. Finally, by an approximation argument, we have that (9) holds for any *<sup>φ</sup>* <sup>∈</sup> *<sup>C</sup>*<sup>∞</sup> *<sup>c</sup>* ((0, *<sup>T</sup>*) <sup>×</sup> <sup>T</sup>3) and we conclude.

#### **4. Approximate Solutions of the Incompressible Navier-Stokes Equations**

In this section, we define the notion of *approximate sequence of solutions* to the Navier-Stokes equations and we prove the convergence to Leray-Hopf weak solutions. We use an approach which is a little different from the one usual used. Our choice, which does not follows the historical path, is motivated by the pedagogical purpose of having a unified treatment for several different methods.

**Definition 2.** *Let <sup>n</sup>* <sup>∈</sup> <sup>N</sup>*. We say that* {*un*}*<sup>n</sup>* <sup>⊂</sup> *<sup>C</sup>*(0, *<sup>T</sup>*; *<sup>L</sup>*2(Ω)) *is an* approximate sequence of solutions *with divergence-free initial datum u<sup>n</sup>* <sup>0</sup> *if*

*1. It holds that*

$$L\{\boldsymbol{u}^{n}\}\_{n} \text{ is a bounded sequence in } L^{\infty}(0,T;H) \cap L^{2}(0,T;V);\tag{10}$$

*2. For any <sup>n</sup>* <sup>∈</sup> <sup>N</sup> *and any <sup>φ</sup>* ∈ V(Ω) *there exists <sup>R</sup><sup>n</sup> <sup>φ</sup>* <sup>∈</sup> *<sup>L</sup>*1(0, *<sup>T</sup>*) *such that for any <sup>χ</sup>* <sup>∈</sup> *C*<sup>∞</sup> *<sup>c</sup>* ([0, *T*))

$$\begin{aligned} \int\_{0}^{T} \left( (\boldsymbol{u}^{n}(t), \boldsymbol{\phi}) \boldsymbol{\chi}(t) + ((\boldsymbol{u}^{n} \cdot \nabla) \boldsymbol{u}^{n}, \boldsymbol{\phi}) \boldsymbol{\chi}(t) + \nu (\nabla \boldsymbol{u}^{n}, \nabla \boldsymbol{\phi}) \boldsymbol{\chi}(t) \right) dt \\ - (\boldsymbol{u}^{n}\_{0}, \boldsymbol{\phi}) \boldsymbol{\chi}(0) = \int\_{0}^{T} \boldsymbol{R}^{n}\_{\boldsymbol{\phi}}(t) \boldsymbol{\chi}(t) \, dt; \end{aligned} \tag{11}$$

*3. It holds*

$$R^{\rm II}\_{\phi} \rightharpoonup 0 \qquad \text{weakly in } L^1(0, T) \text{ as } n \to \infty; \tag{12}$$

*4. For any n* <sup>∈</sup> <sup>N</sup> *and t* <sup>∈</sup> (0, *<sup>T</sup>*) *it holds that*

$$\|\|u^n(t)\|\|\_{2}^{2} + 2\nu \int\_{0}^{t} \|\nabla u^n(s)\|\|\_{2}^{2} \, ds \le \|\|u\_0^n\|\|\_{2}^{2}.\tag{13}$$

Since generally the existence of (smooth) approximating sequences is *rather easy* to be proved, the advantage of this definition is that one has just to check a condition on the data and condition (12) on the remainder (commutator) to show that the approximate solutions converge to a Leray-Hopf weak solution, as is done in the next theorem.

**Theorem 1.** *Let <sup>u</sup>*<sup>0</sup> <sup>∈</sup> *<sup>H</sup> and* {*un*}*<sup>n</sup> be a sequence of approximate solutions with initial data* {*u<sup>n</sup>* 0 }*n such that*

$$
u\_0^n \to \mathfrak{u}\_0 \qquad\text{strongly in } H.\tag{14}$$

*Then, up to a sub-sequence not relabelled, there exists u such that if* Ω *satisfies* (*A*1) *and* (*A*2) *then*

$$
\mu^n \to \mu \qquad\qquad\text{strongly in } L^2(0, T; L^2(\Omega)),\tag{15}
$$

*and if* Ω = R<sup>3</sup>

$$
\mu^n \to \mu \qquad\qquad\text{strongly in } L^2(0, T; L^2\_{loc}(\mathbb{R}^3)).\tag{16}
$$

*Moreover, u is a Leray-Hopf weak solution of* (1) *and* (2)*.*

**Remark 3.** *We stress that because of Remark 1 requiring condition* (13) *(that is a good energy balance already on the approximate functions) is fundamental in order to obtain the energy inequality* (7)*. Moreover, by inspecting the proof below it will be clear that given* {*un*}*<sup>n</sup> satisfying* (1) *and* (2) *in Definition 2 then there exists u satisfying* (1) *and* (2) *in Definition 1 such that the convergences* (15) *and* (16) *hold. This remark will be important in the analysis of the Implicit Euler Scheme in Section 5.3, because the scheme will not fully fit in the framework of Definition 2.*

**Remark 4.** *We point out that Theorem 1 is not limited to the case of vanishing external force, but the result holds also in the presence of an external body force <sup>f</sup>* <sup>∈</sup> *<sup>L</sup>*2(0, *<sup>T</sup>*;(*H*<sup>1</sup> <sup>0</sup> (Ω))- )*. This can be obtained with minor changes in the proof. In this case Definition 2 must also be integrated with an approximating sequence f n, and adding the term*

$$\int\_{0}^{T} (f^n, u^n) \chi(t) \, dt.$$

*to the left-hand side of* (11)*. This is due to the fact that some of the approximation methods in Section <sup>5</sup> require the body force to be smooth and it is enough to require <sup>f</sup> <sup>n</sup>* <sup>→</sup> *<sup>f</sup> in L*2(0, *T*;(*H*<sup>1</sup> <sup>0</sup> (Ω))- )*.*

*We also note that requiring the body force to be in L*2(0, *T*;(*H*<sup>1</sup> <sup>0</sup> (Ω))- ) *(and not only in L*2(0, *T*; *V*- )*) is needed in order to remain inside the space of distributions and to have an associated pressure as in Lemma 2; We refer to Reference [17] for more details on this issue.*

We start with the following straightforward corollary of the classical Arzelà-Ascoli theorem for real functions of a real variable.

**Lemma 3.** *Let E be a separable Banach space and let* E ⊂ *E be a dense subset. Let* {*Fn*}*<sup>n</sup> be a sequence of measurable functions such that Fn* : [0, *T*] → *E*∗*. Assume that*


*Then, Fn* ∈ *Cw*([0, *T*]; *E*∗) *and there exists F* ∈ *Cw*([0, *T*]; *E*∗) *such that, up to a sub-sequence,*

$$F\_n \to F \quad \text{in } \mathbb{C}\_w([0, T]; E^\*).$$

A fundamental step in the proof of existence for nonlinear partial differential equations is the proof of certain compactness which allows to get strong convergence in suitable norms. Observe that the a-priori bounds are useful to get weak or weak-\* convergences, by means of application of the Riesz representation theorem and –more generally– of Banach-Alaoglu-Bourbaki theorem. On the other hand, since *T*(*xn*) - *T*(*x*) for a linear operator *T*, weak convergence allows to consider linear equations, or more precisely, the linear terms in the equations. On the other hand weak convergence is in general not enough to prove that

$$\int\_{0}^{T} (\left(\boldsymbol{u}^{\boldsymbol{n}} \cdot \nabla\right) \boldsymbol{u}^{\boldsymbol{n}}, \boldsymbol{\phi}) \chi(t) \, dt \to \int\_{0}^{T} (\left(\boldsymbol{u} \cdot \nabla\right) \boldsymbol{u}, \boldsymbol{\phi}) \chi(t) \, dt, \qquad \text{as } \boldsymbol{n} \to \infty.$$

Hence, by the a priori estimates we can construct a limit object *u*, but we still have to show that *u* is a weak solution of the limiting problem.

To address this point several results have been used. Leray used Helly's theorem on monotone functions and an ingenious application of Riesz theorem with multiple Cantor diagonal arguments. Hopf used an inequality by Friederichs to handle the Galerkin case. Starting from the work of J.L. Lions [18] it became common to use the approach by the socalled Aubin-Lions lemma, which is borrowed from the general theory of abstract equations and is based on obtaining some estimates on the time derivative (at least in negative space) of the solution. This latter approach is very flexible, but it requires some non-trivial functional analysis preliminaries to estimate the time-derivative, since instead one can use directly some properties coming from the proper definition of the approximation. Note that in the Definition 1 of weak solution there is no mention to the time-derivative. We will show how to obtain compactness in a elementary way, directly from the weak formulation and thus avoiding the use of time derivatives in Bochner spaces. We believe this may be a simpler approach, at least for presentation to students. We also point out that in certain applications to more complex fluid problems as for instance fluids in a moving domain

or non-Newtonian fluids with rheology with time-dependent constitutive law, the proper definition of the time derivative is technically complicated and an approach avoiding the use of this notion becomes particularly welcome.

The next lemma provides a general criterion for strong convergence, which has the advantage to avoid assumptions on the time derivative. Of course, the lemma holds only on bounded domains and therefore we exclude the whole space case, since in the latter one has to work locally. We also stress that the hypothesis are not optimal since we do not prove an *if and only if*, but the hypotheses are easily verifiable for general nonlinear evolution problems.

The lemma below is very similar to the one proved by Landes and Mustonen in Reference [19] and for an application to the Navier-Stokes equations see Landes [20]. For an optimal version (at least in general Hilbert spaces) we refer to Rakotoson and Temam [21].

**Lemma 4.** *Let <sup>U</sup>* <sup>⊂</sup> <sup>R</sup><sup>3</sup> *be any bounded domain or <sup>U</sup>* <sup>=</sup> <sup>T</sup>3*. Let* <sup>1</sup> <sup>&</sup>lt; *<sup>p</sup>* <sup>&</sup>lt; <sup>∞</sup> *and assume that <sup>g</sup>* <sup>∈</sup> *<sup>L</sup>*∞(0, *<sup>T</sup>*; *<sup>L</sup>*1(*U*)) <sup>∩</sup> *<sup>L</sup>p*(0, *<sup>T</sup>*; *<sup>W</sup>*1,*<sup>p</sup>* <sup>0</sup> (*U*)) *and* {*gn*}*<sup>n</sup> is a sequence such that*

> {*gn*}*<sup>n</sup> is bounded in L*∞(0, *<sup>T</sup>*; *<sup>L</sup>*1(*U*)) <sup>∩</sup> *<sup>L</sup>p*(0, *<sup>T</sup>*; *<sup>W</sup>*1,*<sup>p</sup>* <sup>0</sup> (*U*)),

*and gn*(*t*) *<sup>g</sup>*(*t*) *weakly in L*1(*U*) *for a.e. t* <sup>∈</sup> [0, *<sup>T</sup>*]*. Then, it holds that <sup>g</sup><sup>n</sup>* <sup>→</sup> *g in Lp*(0, *<sup>T</sup>*; *<sup>L</sup>p*(*U*)).

**Proof.** We prove the lemma only in the case of *<sup>U</sup>* <sup>⊂</sup> <sup>R</sup><sup>3</sup> being a bounded domain with smooth boundary. First, since *<sup>g</sup>n*(*t*, ·) and *<sup>g</sup>*(*t*, ·) are in *<sup>W</sup>*1,*<sup>p</sup>* <sup>0</sup> (*U*), their extensions to zero off *<sup>U</sup>* are both in *<sup>W</sup>*1,*p*(R3). We denote by *<sup>g</sup>*¯*<sup>n</sup>* and *<sup>g</sup>*¯ these extensions. a.e. *<sup>t</sup>* <sup>∈</sup> (0, *<sup>T</sup>*).

Let *ρε* be a standard spatial mollifier and set *g<sup>n</sup> <sup>ε</sup>* :<sup>=</sup> *ρε* <sup>∗</sup> *<sup>g</sup>*¯*<sup>n</sup>* and *<sup>g</sup><sup>ε</sup>* :<sup>=</sup> *ρε* <sup>∗</sup> *<sup>g</sup>*¯. Next, we have that

$$\begin{split} \left| \mathcal{J}\_{\varepsilon}^{\boldsymbol{u}}(t,\boldsymbol{x}) - \mathcal{J}^{\boldsymbol{u}}(t,\boldsymbol{x}) \right| &\leq \varepsilon \int\_{B\_{1}} \rho(\boldsymbol{y}) \int\_{0}^{1} \left| \nabla \mathcal{J}^{\boldsymbol{u}}(t,\boldsymbol{x} - \boldsymbol{\varepsilon}\boldsymbol{\tau}\boldsymbol{y}) \right| d\boldsymbol{\tau} d\boldsymbol{y} \\ &\leq \varepsilon \left( \int\_{B\_{1}} \rho^{\boldsymbol{p}'}(\boldsymbol{y}) \, d\boldsymbol{y} \right)^{\frac{1}{p'}} \left( \int\_{B\_{1}} \left( \int\_{0}^{1} |\nabla \mathcal{J}^{\boldsymbol{u}}(t,\boldsymbol{x} - \boldsymbol{\varepsilon}\boldsymbol{\tau}\boldsymbol{y})| \, d\boldsymbol{\tau} \right)^{p} \, d\boldsymbol{y} \right)^{\frac{1}{p}} \\ &\leq \varepsilon \Big( \int\_{B\_{1}} \rho^{\boldsymbol{p}'}(\boldsymbol{y}) \, d\boldsymbol{y} \Big)^{\frac{1}{p'}} \left( \int\_{B\_{1}} \int\_{0}^{1} |\nabla \mathcal{J}^{\boldsymbol{u}}(t,\boldsymbol{x} - \boldsymbol{\varepsilon}\boldsymbol{\tau}\boldsymbol{y})|^{p} \, d\boldsymbol{\tau} d\boldsymbol{y} \right)^{\frac{1}{p}}, \end{split}$$

and the same estimate holds also for *g*¯. Therefore, the following estimates hold

$$\begin{aligned} \int\_0^T ||\underline{\boldsymbol{g}}\_{\varepsilon}^{\boldsymbol{\mathsf{u}}} - \underline{\boldsymbol{\mathsf{g}}}^{\boldsymbol{\mathsf{u}}}||\_p^p \, dt &\lesssim \varepsilon^p \int\_0^T ||\nabla \underline{\boldsymbol{g}}^{\boldsymbol{\mathsf{u}}}||\_p^p \, dt, \\ \int\_0^T ||\underline{\boldsymbol{g}}\_{\varepsilon} - \underline{\boldsymbol{\mathsf{g}}}||\_p^p \, dt &\lesssim \varepsilon^p \int\_0^T ||\nabla \underline{\boldsymbol{g}}||\_p^p \, dt, \end{aligned}$$

with bounds depending only on *ρ*.

Next, by triangular inequality we have

$$\begin{split} \int\_{0}^{T} \left\| \boldsymbol{\varrho}^{u} - \boldsymbol{g} \right\|\_{L^{p}(\mathcal{U})}^{p} dt &\leq \int\_{0}^{T} \left\| \boldsymbol{\varrho}^{u} - \boldsymbol{g}\_{\varepsilon}^{u} \right\|\_{L^{p}(\mathcal{U})}^{p} dt + \int\_{0}^{T} \left\| \boldsymbol{g}\_{\varepsilon}^{u} - \boldsymbol{g}\_{\varepsilon} \right\|\_{L^{p}(\mathcal{U})}^{p} dt + \int\_{0}^{T} \left\| \boldsymbol{g}\_{\varepsilon} - \boldsymbol{g} \right\|\_{L^{p}(\mathcal{U})}^{p} dt \\ &\leq \int\_{0}^{T} \left\| \boldsymbol{\varrho}^{u} - \boldsymbol{g}\_{\varepsilon}^{u} \right\|\_{L^{p}(\mathbb{R}^{3})}^{p} dt + \int\_{0}^{T} \left\| \boldsymbol{g}\_{\varepsilon}^{u} - \boldsymbol{g}\_{\varepsilon} \right\|\_{L^{p}(\mathcal{U})}^{p} dt + \int\_{0}^{T} \left\| \boldsymbol{g}\_{\varepsilon} - \boldsymbol{g} \right\|\_{L^{p}(\mathbb{R}^{3})}^{p} dt \\ &\leq \varepsilon^{p} \int\_{0}^{T} \left( \left\| \boldsymbol{\nabla} \boldsymbol{\varrho}^{u} \right\|\_{p}^{p} + \left\| \boldsymbol{\nabla} \boldsymbol{\varrho} \right\|\_{p}^{p} \right) dt + \int\_{0}^{T} \int\_{\mathcal{U}} |\boldsymbol{g}\_{\varepsilon}^{u}(t,x) - \boldsymbol{g}\_{\varepsilon}(t,x)|^{p} \,d\mathbf{x} dt. \end{split} \tag{17}$$

The first term from the right-hand side can be made arbitrarily small by choosing *ε* small enough.

To conclude, we first note that by definition of convolution, there exists *C* = *C*(*ε*, *U*) such that

$$|\mathcal{g}\_{\varepsilon}(t, \mathbf{x})| + |\mathcal{g}\_{\varepsilon}^{n}(t, \mathbf{x})| \le \mathbf{C}.$$

Next, since clearly it holds that for the extended functions it holds

$$
\bar{\mathfrak{g}}^{\mathfrak{n}}(t) \rightharpoonup \bar{\mathfrak{g}}(t) \quad \text{weakly in } L^1(\mathbb{R}^3),
$$

we also have that

$$g\_{\varepsilon}^{\mathfrak{n}}(t, \mathfrak{x}) \to g\_{\varepsilon}(t, \mathfrak{x}), \quad \text{a.e. in } t \in (0, T) \text{ and for all } \mathfrak{x} \in \mathbb{R}^3.$$

This follows by fixing > 0, a time *t* such that *g*¯*n*(*t*) *<sup>g</sup>*¯(*t*), *<sup>x</sup>* <sup>∈</sup> <sup>R</sup>3, and noticing that

$$
\left(\underline{\mathcal{g}}\_{\varepsilon}^{\eta}(t,\mathtt{x}) - \overline{\mathcal{g}}\_{\varepsilon}(t,\mathtt{x})\right) = \int\_{\mathbb{R}^{3}} (\underline{\mathcal{g}}^{\eta}(t,\mathtt{y}) - \underline{\mathcal{g}}(t,\mathtt{y})) \rho\_{\varepsilon}(\mathtt{x} - \boldsymbol{y}) \, d\boldsymbol{y} \to 0,
$$

as *n* → ∞.

This shows that, for any fixed *ε* > 0, the last term in last inequality in (17) goes to zero as *n* → ∞, by using Dominated Convergence Theorem. The proof is concluded since we showed that *<sup>g</sup><sup>n</sup>* <sup>−</sup> *<sup>g</sup> <sup>L</sup>p*(0,*T*;*Lp*(*U*)) can be made arbitrarily small.

The following theorem is the main result of this section.

**Proof of the Theorem 1.** Let {*un*}*<sup>n</sup>* be a sequence of approximate solutions. By condition (10) of Definition 2 we can infer that up to a sub-sequence (not relabelled) there exists *<sup>u</sup>* <sup>∈</sup> *<sup>L</sup>*∞(0, *<sup>T</sup>*; *<sup>H</sup>*) <sup>∩</sup> *<sup>L</sup>*2(0, *<sup>T</sup>*; *<sup>V</sup>*) such that

$$
u^n \stackrel{\ast}{\rightharpoonup} u \qquad \text{weakly-\* in } L^\infty(0, T; L^2(\Omega)), \tag{18}$$

$$
u^n \rightharpoonup u \qquad \text{weakly in } L^2(0, T; L^2(\Omega)), \tag{19}$$

$$
\nabla \mu^n \rightharpoonup \nabla \mu \qquad \text{weakly in } L^2(0, T; L^2(\Omega)). \tag{20}
$$

For *<sup>k</sup>* <sup>∈</sup> <sup>N</sup> and ¯*<sup>t</sup>* <sup>∈</sup> (0, *<sup>T</sup>*), we consider the following function

$$\chi\_k^{\mathbb{F}}(t) = \begin{cases} 1, & t \in [0, \overline{t}) \\ k(\overline{t} - t) + 1, & t \in [\overline{t}, \overline{t} + \frac{1}{k}) \\ 0, & t \in [\overline{t} + \frac{1}{k}, T), \end{cases} \tag{21}$$

Let *<sup>φ</sup>* <sup>∈</sup> *<sup>C</sup>*<sup>∞</sup> *<sup>c</sup>* (Ω) with div *<sup>φ</sup>* <sup>=</sup> 0 and *<sup>s</sup>*, *<sup>t</sup>* <sup>∈</sup> (0, *<sup>T</sup>*). By using the function *<sup>χ</sup>*¯*<sup>t</sup> <sup>k</sup>*(*t*) with first with ¯*<sup>t</sup>* <sup>=</sup> *<sup>t</sup>* and then with ¯*<sup>t</sup>* <sup>=</sup> *<sup>s</sup>*, together with the fact that *<sup>u</sup><sup>n</sup>* <sup>∈</sup> *<sup>C</sup>*(0, *<sup>T</sup>*; *<sup>L</sup>*2(Ω)) we can infer that

$$\left( \left( u^n(t), \phi \right) - \left( u^n(s), \phi \right) + \int\_s^t \left( \left( \left( u^n(\tau) \cdot \nabla \right) u^n(\tau), \phi \right) + \nu \left( \nabla u^n(\tau), \nabla \phi \right) + R^n\_{\phi}(\tau) \right) d\tau = 0. \right)$$

Next, for *<sup>φ</sup>* <sup>∈</sup> *<sup>C</sup>*<sup>∞</sup> *<sup>c</sup>* (Ω), let *<sup>F</sup>n*(*t*) := ((*un*(*t*) · ∇) *<sup>u</sup>n*(*t*), *<sup>φ</sup>*) + *<sup>ν</sup>*(∇*un*(*t*), <sup>∇</sup>*φ*) + *<sup>R</sup><sup>n</sup> <sup>φ</sup>*(*t*). Then, condition (10), the Gagliardo-Nirenberg-Sobolev inequality (4), and the hypothesis on *R<sup>n</sup> <sup>φ</sup>* in (11) imply that the family {*Fn*}*<sup>n</sup>* is equi-integrable and then the function *<sup>t</sup>* → (*un*(*t*), *<sup>φ</sup>*) is equi-continuous. Since <sup>V</sup>(Ω) is dense in *<sup>H</sup>* and *<sup>H</sup>* is reflexive, we can conclude by using Lemma 3 that

$$
\mu^n \to \mu \qquad \text{in } \mathbb{C}\_w([0, T]; H). \tag{22}
$$

By using (22) and (20) we can prove that *u* satisfies the energy inequality (7). Indeed, for any *t* ∈ (0, *T*) we have that

$$\begin{aligned} \|\boldsymbol{u}(t)\|\|\_{2}^{2} &\leq \liminf\_{n\to\infty} \|\boldsymbol{u}^{n}(t)\|\|\_{2'}^{2} \\ \int\_{0}^{t} \|\nabla\boldsymbol{u}(s)\|\|\_{2}^{2} \, ds &\leq \liminf\_{n\to\infty} \int\_{0}^{t} \|\nabla\boldsymbol{u}^{n}(s)\|\|\_{2}^{2} \, ds. \end{aligned}$$

Then,

$$\begin{split} \|\|u(t)\|\|\_{2}^{2} + 2\nu \int\_{0}^{t} \|\nabla u(s)\|\|\_{2}^{2} ds \leq & \mathop{\rm lim}\limits\_{n\to\infty} \|\|u^{n}(t)\|\|\_{2}^{2} + \mathop{\rm lim}\limits\_{n\to\infty} \|\operatorname{2}\nu \int\_{0}^{t} \|\nabla u^{n}(s)\|\|\_{2}^{2} ds \\ \leq & \mathop{\rm lim}\limits\_{n\to\infty} \Big( \|\|u^{n}(t)\|\|\_{2}^{2} + 2\nu \int\_{0}^{t} \|\nabla u^{n}(s)\|\|\_{2}^{2} ds \Big) \\ \leq & \mathop{\rm lim}\limits\_{n\to\infty} \|u\_{0}^{n}\|\|\_{2}^{2} = \|u\_{0}\|\|\_{2}^{2} .\end{split}$$

where we have used (14). In order to conclude it remains only to prove (15) and (16). If Ω is the flat torus or a bounded domain, then (15) follows directly by Lemma 4. If Ω = R<sup>3</sup> we need a localization argument. We first note that by the Helmholtz decomposition (22) holds also in *Cw*([0, *<sup>T</sup>*]; *<sup>L</sup>*2(R3)). Next, for *<sup>k</sup>* <sup>∈</sup> <sup>N</sup> let *<sup>ψ</sup>* <sup>∈</sup> *<sup>C</sup>*<sup>∞</sup> *<sup>c</sup>* (*Bk*<sup>+</sup>1(0)) such that *ψ* = 1 on *Bk* and define *<sup>g</sup><sup>n</sup>* :<sup>=</sup> *<sup>u</sup><sup>n</sup> <sup>ψ</sup>*. The sequence {*gn*}*<sup>n</sup>* satisfies the hypothesis of Lemma 4, and therefore, after a diagonal argument, it follows that there exists a sub-sequence not relabelled such that

$$\mathfrak{g}'' \to \mathfrak{g} = \mathfrak{u}\,\psi \qquad \text{strongly in } L^2(0, T; L^2(B\_{k+1})).\tag{23}$$

Then, condition (23) easily implies (16).

#### **5. Approximation Methods**

After the general result of the previous section, we are now going to show that a general class of methods used to construct weak solutions will fit the in the framework of Theorem 1, as described in Section 4.

#### *5.1. Leray Approximation Scheme*

We start describing the original scheme introduced by Leray in Reference [2] (even if we use a completely different compactness argument to show the convergence of approximations). In this case we consider <sup>Ω</sup> <sup>=</sup> <sup>R</sup>3. We fix a sequence {*εn*} of positive numbers going to zero and let *ρε<sup>n</sup>* be a standard mollifier (only) in the space variables. For *<sup>v</sup>* : (0, *<sup>T</sup>*) <sup>×</sup> <sup>R</sup><sup>3</sup> → <sup>R</sup><sup>3</sup> we set <sup>Ψ</sup>*n*(*v*) :<sup>=</sup> *ρε<sup>n</sup>* <sup>∗</sup> *<sup>v</sup>*, where the convolution is only in the space variables. Let *<sup>u</sup>*<sup>0</sup> <sup>∈</sup> *<sup>H</sup>* and let *<sup>n</sup>* <sup>∈</sup> <sup>N</sup>. Define *<sup>u</sup><sup>n</sup>* <sup>0</sup> = <sup>Ψ</sup>*n*(*u*0) and *<sup>u</sup><sup>n</sup>* as the solution of the following Cauchy problem:

$$\begin{cases} \partial\_t u^n - \nu \Delta u^n + \left( \mathbb{V}\_n(u^n) \cdot \nabla \right) u^n + \nabla p^n = 0 & \text{in } (0, T) \times \mathbb{R}^3 \\\\ \text{div} \, u^n = 0 & \text{in } (0, T) \times \mathbb{R}^3 \\\\ u^n|\_{t=0} = u\_0^n & \text{on } \{t = 0\} \times \mathbb{R}^3 \end{cases} \tag{24}$$

We want to prove that for any *<sup>n</sup>* <sup>∈</sup> <sup>N</sup> the function *<sup>u</sup><sup>n</sup>* exists, is smooth, and {*un*}*<sup>n</sup>* is an approximate sequence of solutions in the sense of Definition 2.

**Theorem 2.** *Let u*<sup>0</sup> ∈ *H. Then, it holds that*


$$u^n \to u \qquad \text{strongly in } L^2(0, T; L^2\_{loc}(\mathbb{R}^3)).$$

**Proof.** Let us prove (1). The proof is very classical so we only sketch it. By using a fixed point argument we can prove that there exists a time *<sup>T</sup>*<sup>1</sup> <sup>=</sup> *<sup>T</sup>*1( *<sup>u</sup><sup>n</sup>* <sup>0</sup> *<sup>H</sup>*<sup>3</sup> ) > 0 such that there exists a unique *<sup>u</sup><sup>n</sup>* <sup>∈</sup> *<sup>C</sup>*([0, *<sup>T</sup>*∗); *<sup>H</sup>*3(R3)) solution of (24) for *<sup>T</sup>*<sup>∗</sup> <sup>≥</sup> *<sup>T</sup>*1. Let us suppose that *<sup>T</sup>*<sup>∗</sup> is the maximal time of existence of *<sup>u</sup><sup>n</sup>* and, if *<sup>T</sup>*<sup>∗</sup> <sup>&</sup>lt; *<sup>T</sup>* then lim*t*→*T*∗− *<sup>u</sup>n*(*s*) *<sup>H</sup>*<sup>3</sup> <sup>=</sup> <sup>∞</sup>.

To obtain a global solution we exploit a standard energy estimates argument. Indeed, we first note that by multiplying (24) by *u<sup>n</sup>* and integrating by parts we get

$$||\boldsymbol{u}^{n}(t)||\_{2}^{2} + 2\nu \int\_{0}^{t} ||\nabla \boldsymbol{u}^{n}(s)||\_{2}^{2} ds = ||\boldsymbol{u}\_{0}^{n}||\_{2\prime}^{2} \tag{25}$$

with an equality which is valid for all *t* < *T*∗. Note that this is exactly the same calculation we have done formally to obtain the energy inequality (3) in the introduction. In particular, from (25) we obtain

$$\sup\_{t \in (0, T)} \|u^\eta(t)\|\_2 \le \|u\_0^\eta\|\_2. \tag{26}$$

Next, by using that *H*3(R3) is an algebra and by using the standard properties of mollifiers, it is easy to prove that

$$\begin{split} \frac{d}{dt} \|\boldsymbol{u}^{n}(t)\|\_{H^{3}}^{2} + \nu \|\nabla\boldsymbol{u}^{n}(t)\|\_{H^{3}}^{2} &\lesssim \left| \left( (\boldsymbol{\Psi}\_{n}(\boldsymbol{u}^{n}(t)) \cdot \nabla) \boldsymbol{u}^{n}(t), \boldsymbol{u}^{n}(t) \right) \right|\_{H^{3}} \| \\ &\lesssim \|\boldsymbol{\Psi}\_{n}(\boldsymbol{u}^{n}(t))\|\_{H^{3}} \|\nabla\boldsymbol{u}^{n}(t)\|\_{H^{3}} \|\boldsymbol{u}^{n}(t)\|\_{H^{3}} \\ &\lesssim \|\boldsymbol{\Psi}\_{n}(\boldsymbol{u}^{n}(t))\|\_{H^{3}}^{2} \|\boldsymbol{u}^{n}(t)\|\_{H^{3}}^{2} + \frac{\nu}{2} \|\nabla\boldsymbol{u}^{n}(t)\|\_{H^{3}}^{2} \\ &\lesssim \frac{1}{\epsilon\_{n}^{6}} \|\boldsymbol{u}^{n}(t)\|\_{2}^{2} \|\boldsymbol{u}^{n}(t)\|\_{H^{3}}^{2} + \frac{\nu}{2} \|\nabla\boldsymbol{u}^{n}(t)\|\_{H^{3}}^{2} \\ &\lesssim \frac{1}{\epsilon\_{n}^{6}} \|\boldsymbol{u}\_{0}^{n}\|\_{2}^{2} \|\boldsymbol{u}^{n}(t)\|\_{H^{3}}^{2} + \frac{\nu}{2} \|\nabla\boldsymbol{u}^{n}(t)\|\_{H^{3}}^{2}. \end{split}$$

where in the last inequality we have used (26). Therefore, we have that

$$\frac{d}{dt} \|u^n(t)\|\_{H^3}^2 \le C\_n \|u\_0^n\|\_2^2 \|u^n(t)\|\_{H^3}^2$$

and, by using the Gronwall Lemma, we conclude that necessarily *T*∗ = *T* (this argument shows that in fact *<sup>u</sup><sup>n</sup>* is defined for all *<sup>t</sup>* <sup>&</sup>gt; 0, for any fixed *<sup>n</sup>* <sup>∈</sup> <sup>N</sup>).

Next, to show (2) it is enough to prove that {*un*}*<sup>n</sup>* satisfies the conditions in Definition 2. Clearly, from(1) we have that {*un*}*<sup>n</sup>* <sup>⊂</sup> *<sup>C</sup>*([0, *<sup>T</sup>*]; *<sup>L</sup>*2(R3)). By using the standard property of mollifiers <sup>Ψ</sup>*n*(*u*0) <sup>2</sup> <sup>≤</sup>*u*<sup>0</sup> 2, and from the energy estimate (25) we get that {*un*}*<sup>n</sup>* is bounded uniformly in *<sup>L</sup>*∞(0, *<sup>T</sup>*; *<sup>H</sup>*) <sup>∩</sup> *<sup>L</sup>*2(0, *<sup>T</sup>*; *<sup>H</sup>*). Moreover, (25) is exactly (13) and then it remains only to verify (11). For *<sup>φ</sup>* <sup>∈</sup> *<sup>C</sup>*<sup>∞</sup> *<sup>c</sup>* (Ω) and *<sup>t</sup>* <sup>∈</sup> (0, *<sup>T</sup>*) the function *<sup>t</sup>* → *<sup>R</sup><sup>n</sup> <sup>φ</sup>*(*t*) is defined as follows

$$\mathcal{R}^n\_{\phi}(t) := \left( \left( \left[ \Psi\_n(\boldsymbol{u}^n(t)) - \boldsymbol{u}^n(t) \right) \right] \cdot \nabla \right) \boldsymbol{u}^n(t), \phi \right).$$

With this choice of *R<sup>n</sup> <sup>φ</sup>*, the equation (11) is satisfied and it remains only to prove that convergence stated in condition (12) of Definition 2. First, note that by Hölder inequality

$$|R^n\_\Phi(t)| \le ||\nabla u^n(t)||\_2 ||\phi||\_\infty ||\Psi\_n(u^n(t)) - u^n(t)||\_2. \tag{27}$$

Then, for any fixed (*t*, *<sup>x</sup>*) <sup>∈</sup> (0, *<sup>T</sup>*) <sup>×</sup> <sup>R</sup>3, by a direct calculation (using again the properties of mollifiers) we have

$$\begin{split} \left| \Psi\_{n} (\boldsymbol{u}^{n} (t, \boldsymbol{x})) - \boldsymbol{u}^{n} (t, \boldsymbol{x}) \right| &\leq \varepsilon\_{n} \int\_{B\_{1}} \rho(\boldsymbol{y}) \int\_{0}^{1} \left| \nabla \boldsymbol{u}^{n} (t, \boldsymbol{x} - \boldsymbol{\varepsilon} \boldsymbol{\tau} \boldsymbol{y}) \right| d\boldsymbol{\tau} d\boldsymbol{y} \\ &\leq \varepsilon\_{n} \left( \int\_{B\_{1}} \rho^{2}(\boldsymbol{y}) \, d\boldsymbol{y} \right)^{\frac{1}{2}} \left( \int\_{B\_{1}} \left( \int\_{0}^{1} |\nabla \boldsymbol{u}^{n} (t, \boldsymbol{x} - \boldsymbol{\varepsilon} \boldsymbol{\tau} \boldsymbol{y})| \, d\boldsymbol{\tau} \right)^{2} \, d\boldsymbol{y} \right)^{\frac{1}{2}} \\ &\leq \varepsilon\_{n} \left( \int\_{B\_{1}} \rho^{2}(\boldsymbol{y}) \, d\boldsymbol{y} \right)^{\frac{1}{2}} \left( \int\_{B\_{1}} \int\_{0}^{1} |\nabla \boldsymbol{u}^{n} (t, \boldsymbol{x} - \boldsymbol{\varepsilon} \boldsymbol{\tau} \boldsymbol{y})|^{2} \, d\boldsymbol{\tau} d\boldsymbol{y} \right)^{\frac{1}{2}}. \end{split}$$

Then, by a further integration

$$\int\_{\mathbb{R}^3} |\Psi\_n(\mu^n(t, \mathbf{x})) - \mu^n(t, \mathbf{x})|^2 \, d\mathbf{x} \le \varepsilon\_n^2 \int\_{\mathbb{R}^2} |\nabla \mu^n(t, \mathbf{x})|^2 \, d\mathbf{x} \, d\mathbf{x}$$

and going back to (27) we have that

$$\left|R^{\boldsymbol{\mathfrak{u}}}\_{\boldsymbol{\Phi}}(t)\right| \leq \varepsilon\_{\boldsymbol{\mathfrak{u}}} \left||\nabla \boldsymbol{\mathfrak{u}}^{\boldsymbol{\mathfrak{u}}}(t)\right||\_{2}^{2} \left\|\boldsymbol{\mathfrak{g}}\right\|\_{\infty}.$$

Since the sequence {∇*un*}*<sup>n</sup>* is bounded uniformly with respect to *<sup>n</sup>* in *<sup>L</sup>*2(0, *<sup>T</sup>*; *<sup>L</sup>*2(R3)), we have that *R<sup>n</sup> <sup>φ</sup>* <sup>→</sup> 0 in *<sup>L</sup>*1(0, *<sup>T</sup>*).

#### *5.2. Faedo-Galerkin Method*

The next scheme we consider is the Faedo-Galerkin method. The variant we present is close to the one considered by Hopf and is at the basis of several computational methods, which are used also in fields different from fluid dynamics. In particular, we will see that the unified treatment is possible under the assumption of having a basis which is orthogonal in both *L*<sup>2</sup> and *H*1, as is the case of the spectral basis made by eigenfunctions of the Stokes operator. Observe that in the space-periodic case this basis is explicitly constructed by considering complex exponentials, while in the case of a smooth bounded domain, the existence is obtained via the standard theory of compact operators, showing existence of countable non-decreasing positive {*λj*} and smooth {*ψj*} such that it holds for all *<sup>j</sup>* <sup>∈</sup> <sup>N</sup>

$$\begin{aligned} -\Delta \psi\_j + \nabla \pi\_j &= \lambda\_j \psi\_j \quad &\text{in } \Omega, \\ \text{div } \psi\_j &= 0 \quad &\text{in } \Omega, \\ \psi\_j &= 0 \quad &\text{on } \partial \Omega. \end{aligned}$$

We consider <sup>Ω</sup> <sup>⊂</sup> <sup>R</sup><sup>3</sup> with smooth boundary *<sup>∂</sup>*<sup>Ω</sup> or the three-dimensional flat torus. Let be given an orthonormal basis {*ψm*}*m*∈<sup>N</sup> of *<sup>H</sup>*, such that *<sup>ψ</sup><sup>m</sup>* ∈ V(Ω). The Faedo-Galerkin method is based on the construction of approximate solutions of the type

$$\mu^n(t, x) = \sum\_{j=1}^n c\_j^n(t)\psi\_j(x) \qquad n \in \mathbb{N},\tag{28}$$

which solve the Navier-Stokes equations projected equations over the finite dimensional space *Vn* <sup>=</sup> Span(*ψ*1, ... , *<sup>ψ</sup>n*) <sup>⊂</sup> *<sup>V</sup>*. This means that for *<sup>n</sup>* <sup>∈</sup> <sup>N</sup>, the approximate problem to be solved is given by

$$\begin{cases} \frac{d}{dt}(\boldsymbol{u}^{\boldsymbol{n}},\boldsymbol{\psi}\_{\boldsymbol{m}}) + \nu(\nabla\boldsymbol{u}^{\boldsymbol{n}},\nabla\boldsymbol{\psi}\_{\boldsymbol{m}}) + (\left(\boldsymbol{u}^{\boldsymbol{n}}\cdot\nabla\right)\boldsymbol{u}^{\boldsymbol{n}},\boldsymbol{\psi}\_{\boldsymbol{m}}) = \boldsymbol{0} & \boldsymbol{t} \in (0,T),\\ (\boldsymbol{u}^{\boldsymbol{n}}(0),\boldsymbol{\Psi}\_{\boldsymbol{m}}) = (\boldsymbol{u}\_{0\prime}\boldsymbol{\Psi}\_{\boldsymbol{m}}) & \boldsymbol{t} = \boldsymbol{0}, \end{cases} \tag{29}$$

for *m* = 1, ... , *n*, which is a Cauchy problem for a system of *n* ODE's in the coefficients {*cn <sup>j</sup>* (*t*)}*<sup>n</sup> <sup>j</sup>*=1. Let *<sup>P</sup><sup>n</sup>* be projection operator from *<sup>H</sup>* into *Vn*:

$$P^n \colon f \in H \mapsto P^n f := \sum\_{m=1}^n \left( f, \psi\_m \right) \psi\_m \dots$$

Then, the ODE's (29) reduce to the following system of PDE's:

$$\begin{cases} \partial\_t u^n + P^\mu \left( \left( u^n \cdot \nabla \right) u^n \right) - \nu \Delta u^n = 0 & \text{in } (0, T) \times \Omega, \\ \qquad u^n|\_{t=0} = P^n u^0 & \text{in } \Omega. \end{cases} \tag{30}$$

In the next theorem we prove that *<sup>u</sup><sup>n</sup>* is smooth and exists on (0, *<sup>T</sup>*), and that {*un*}*<sup>n</sup>* is an approximate sequence of solutions.

**Theorem 3.** *Let u*<sup>0</sup> ∈ *H. Then, it holds that*


$$
\mu^n \to \mathfrak{u} \qquad \text{strongly in } L^2(0, T; L^2(\Omega)).
$$

**Proof.** We prove (1). By the theory of ordinary differential equations one easily obtains that there exists a unique solution *c<sup>n</sup> <sup>j</sup>* (*t*) <sup>∈</sup> *<sup>C</sup>*1(0, *<sup>T</sup>n*), for some 0 <sup>&</sup>lt; *<sup>T</sup><sup>n</sup>* <sup>≤</sup> *<sup>T</sup>*, being (29) <sup>a</sup> nonlinear (quadratic) system in the coefficients *c<sup>n</sup> <sup>j</sup>* (*t*). Moreover, *<sup>u</sup><sup>n</sup>* is defined through (28) and satisfies (30). Then, by multiplying (30) by *u<sup>n</sup>* and integrating by parts we get

$$\|\|u^n(t)\|\|\_{2}^{2} + 2\nu \int\_{0}^{t} \|\nabla u^n(s)\|\_{2}^{2} ds = \|u\_0^n\|\_{2}^{2} \le \|u\_0\|\_{2\prime}^{2} \tag{31}$$

where we have used that *<sup>u</sup><sup>n</sup>* <sup>0</sup> <sup>2</sup> <sup>=</sup> *<sup>P</sup>nu*<sup>0</sup> <sup>2</sup> <sup>≤</sup>*u*<sup>0</sup> 2. Therefore, for any *<sup>n</sup>* <sup>∈</sup> <sup>N</sup> we have that

$$\sum\_{j=1}^{n} |c\_j^n(t)|^2 = ||\mu^n(t)||\_2^2 \le ||\mu\_0||\_{2'}^2$$

which easily implies that necessarily *T<sup>n</sup>* = *T*.

To prove (2) we show that {*un*}*<sup>n</sup>* satisfy the conditions in Definition 2. Clearly, the sequence {*un*}*<sup>n</sup>* is in *<sup>C</sup>*([0, *<sup>T</sup>*]; *<sup>L</sup>*2(Ω)) and by (31) it verifies the condition (1) and the energy inequality (13). To check that (11) is verified, let *φ* ∈ V(Ω) and, for any *t* ∈ (0, *T*), define

$$R^n\_{\phi}(t) := \left(P^n(\left(\boldsymbol{u}^n(t) \cdot \nabla\right)\boldsymbol{u}^n(t)) - \left(\boldsymbol{u}^n(t) \cdot \nabla\right)\boldsymbol{u}^n(t), \phi\right).$$

Note that we have that

$$\begin{aligned} \|R^n\_{\Phi}(t)\| \lesssim & \|\mu^n(t)\| \|\,\_3\|\|\nabla\mu^n(t)\|\|\_{2}\|\|\Phi - P^n\Phi\|\|\_{\Phi} \\ \lesssim & \|\|\mu^n(t)\|\|\_{2}^{\frac{1}{2}}\|\|\mu^n(t)\|\|\_{H^1}^{\frac{3}{2}}\|\|\Phi - P^n\Phi\|\|\_{H^1\nu} \end{aligned}$$

where we have used the Gagliardo-Nirenberg-Sobolev inequality (4) and that *P<sup>n</sup>* is a projection in both *H* and *V*, since in this case it holds

$$\left\| \left| \nabla \left[ f - \sum\_{m=1}^{n} (f, \psi\_{m}) \psi\_{m} \right] \right| \right\|^{2} = \sum\_{m=n+1}^{\infty} \lambda\_{m} |(f, \psi\_{m})|^{2} \left\| \psi\_{m} \right\|^{2} \stackrel{n \to +\infty}{\longrightarrow} 0,$$

for all *<sup>f</sup>* <sup>∈</sup> *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω). Then, by using Hölder inequality, and Gagliardo-Nirenberg Sobolev inequality and taking into account that *T* < ∞ we have that *R<sup>n</sup> <sup>φ</sup>* <sup>→</sup> 0 in *<sup>L</sup>*1(0, *<sup>T</sup>*).

As already specified if <sup>Ω</sup> <sup>=</sup> <sup>T</sup>3, then one can take {*ψm*}*m*∈<sup>N</sup> to be the Fourier basis. Then, the Faedo-Galerkin method consists in finding the approximated sequence of type (28) solving the Navier-Stokes equations projected over the first *n* Fourier modes. On the other hand, in the case Ω = R<sup>3</sup> one possible choice is to use the method of invading domains, that is to consider the problem in the ball *B*(0, *R*) with zero boundary conditions on *∂B*(0, *R*) and to construct a solution *uR* by the Galerkin method. It turns out that the energy estimate (3) is valid for *uR*, providing uniform estimates (on *uR* which is considered as a function over the whole space, after extension by zero off of Ω); this allows to pass to the limit as *R* → +∞, more or less in the same way as before.

#### *5.3. Implicit Euler Scheme*

The scheme we consider in the present subsection deals with the time-discretization and represents a first step also in the numerical analysis of the Navier-Stokes equations.

We consider the case of Ω being a bounded domain satisfying the hypothesis (*A*3). Let *<sup>n</sup>* <sup>∈</sup> <sup>N</sup> and define the time-step *<sup>κ</sup><sup>n</sup>* :<sup>=</sup> *<sup>T</sup>*/*<sup>n</sup>* and the net *<sup>I</sup><sup>M</sup>* <sup>=</sup> {*tm*}*<sup>n</sup> <sup>m</sup>*=0, such that *ti* − *ti*−<sup>1</sup> = *<sup>κ</sup>n*, for any *<sup>i</sup>* = 1, ..., *<sup>n</sup>*.

Moreover, given *u*<sup>0</sup> ∈ *H*, consider a sequence (In the space periodic setting this can be obtained simply with a mollification with kernel *ρ*, with = <sup>√</sup><sup>1</sup> *<sup>n</sup>* .) of initial data {*u<sup>n</sup>* <sup>0</sup> }*<sup>n</sup>* ⊂ *V* such that

$$\|\|\nabla u\_0^n\|\|\_2 \lesssim \sqrt{n} \|\|u\_0^n\|\|\_{2'} \qquad \text{and} \qquad u\_0^n \to u\_0 \qquad \text{strongly in } H. \tag{32}$$

For *<sup>m</sup>* ∈ {1, ..., *<sup>n</sup>*}, given *<sup>u</sup>*˜*m*−<sup>1</sup> *<sup>n</sup>* the iterate *<sup>u</sup>*\$*<sup>m</sup> <sup>n</sup>* is obtained by solving the boundary value problem

$$\begin{cases} \frac{\hat{\boldsymbol{u}}\_{n}^{m} - \hat{\boldsymbol{u}}\_{n}^{m-1}}{\kappa\_{\text{H}}} - \nu \Delta \hat{\boldsymbol{u}}\_{n}^{m} + \left( \hat{\boldsymbol{u}}\_{n}^{m} \cdot \nabla \right) \hat{\boldsymbol{u}}\_{n}^{m} + \nabla \hat{p}\_{n}^{m} = \boldsymbol{0} & \text{in } \Omega, \\\\ \nabla \cdot \hat{\boldsymbol{u}}\_{n}^{m} = 0 & \text{in } \Omega, \\\\ \hat{\boldsymbol{u}}\_{n}^{m} = 0 & \text{on } \partial \Omega, \end{cases}$$

with *<sup>u</sup>*\$<sup>0</sup> *<sup>n</sup>* = *u<sup>n</sup>* <sup>0</sup> . (In the case of a non-zero force one has to set *f* \$*m <sup>n</sup>* = *κ*−<sup>1</sup> *<sup>n</sup>* ( *tm tm*−<sup>1</sup> *<sup>f</sup>*(*t*) *dt* in the right-hand side of the momentum equation which defines *<sup>u</sup>*\$*<sup>m</sup> <sup>n</sup>* .)

For any fixed *<sup>n</sup>* <sup>∈</sup> <sup>N</sup>, we define the following sequences of functions defined on [0, *<sup>T</sup>*] with values in *V* and in *L*2(Ω):

$$\begin{split} u^{n}(t) &= \sum\_{m=1}^{n} \chi\_{[t\_{m-1}t\_{m}]}(t) \left( \hat{u}\_{n}^{m-1} + \frac{(t - t\_{m-1})}{\kappa\_{n}} (\hat{u}\_{n}^{m} - \hat{u}\_{n}^{m-1}) \right), \quad u^{n}(t\_{n}) = \hat{u}\_{n}^{n}, \\ v^{n}(t) &= \sum\_{m=1}^{n} \chi\_{(t\_{m-1}t\_{m}]}(t) \hat{u}\_{n}^{m}, \quad v^{n}(t\_{0}) = u\_{0}^{n}. \\ p^{n}(t) &= \sum\_{m=1}^{n} \chi\_{(t\_{m-1}t\_{m})}(t) \hat{p}\_{n}^{m}. \end{split} \tag{33}$$

We are now ready to prove the following theorem, which is referred in literature as an "alternate proof" by semi-discretization, see Reference [11].

#### **Theorem 4.** *Let u*<sup>0</sup> ∈ *H. Then, it holds that*

*1. For any fixed <sup>n</sup>* <sup>∈</sup> <sup>N</sup>*, there exist* {*u*\$*<sup>m</sup> <sup>n</sup>* }*<sup>n</sup> <sup>m</sup>*=<sup>1</sup> <sup>⊂</sup> *<sup>H</sup>*<sup>1</sup> <sup>0</sup> (Ω) *such that for any m* = 1, ..., *n and any ψ* ∈ *V*

$$\left(\hat{\boldsymbol{u}}\_{n}^{m},\boldsymbol{\psi}\right) - \left(\hat{\boldsymbol{u}}\_{n}^{m-1},\boldsymbol{\psi}\right) + \kappa\_{n}\boldsymbol{\nu}\left(\nabla\hat{\boldsymbol{u}}\_{n}^{m},\nabla\boldsymbol{\psi}\right) + \kappa\_{n}\left(\left(\hat{\boldsymbol{u}}\_{n}^{m}\cdot\nabla\right)\hat{\boldsymbol{u}}\_{n}^{m},\boldsymbol{\psi}\right) = 0;\tag{34}$$

*2. There exists a Leray-Hopf weak solution u such that the sequence* {*un*}*<sup>n</sup> and* {*vn*}*n, defined in* (33)*, satisfy*

$$\begin{aligned} \boldsymbol{u}^n &\to \boldsymbol{u} &\text{strongly in } L^2(0, T; L^2(\Omega)),\\ \boldsymbol{u}^n - \boldsymbol{v}^n &\to 0 &\text{strongly in } L^2(0, T; L^2(\Omega)). \end{aligned}$$

**Proof.** For the proof of *(1)* we refer to Reference [11]. The idea is the following: For any fixed *<sup>n</sup>* <sup>∈</sup> <sup>N</sup> and any *<sup>m</sup>* <sup>=</sup> 1, ..., *<sup>n</sup>*, the existence of *<sup>u</sup>*\$*<sup>m</sup> <sup>n</sup>* ∈ *V* solution of (34) is obtained by applying the Brouwer fixed point theorem to the following modified version of the steady Navier-Stokes equations, where the given iterate *<sup>u</sup>*\$*m*−<sup>1</sup> *<sup>n</sup>* is considered an external force:

$$\begin{aligned} \frac{\hat{u}\_{\boldsymbol{n}}^{\rm m}}{\kappa\_{\boldsymbol{n}}} - \nu \Delta \hat{u}\_{\boldsymbol{n}}^{\rm m} + \left( \hat{u}\_{\boldsymbol{n}}^{\rm m} \cdot \nabla \right) \hat{u}\_{\boldsymbol{n}}^{\rm m} + \nabla \hat{p}\_{\boldsymbol{n}}^{\rm m} - \frac{\hat{u}\_{\boldsymbol{n}}^{\rm m-1}}{\kappa\_{\boldsymbol{n}}} &= 0 \quad \text{in } \Omega, \\ \nabla \cdot \hat{u}\_{\boldsymbol{n}}^{\rm m} &= 0 \quad \text{in } \Omega, \\ \hat{u}\_{\boldsymbol{n}} &= 0 \quad \text{in } \partial \Omega. \end{aligned}$$

In particular, by the definitions (33), we have that {*un*}*<sup>n</sup>* <sup>⊂</sup> *<sup>C</sup>*(0, *<sup>T</sup>*; *<sup>L</sup>*2(Ω)) and {*vn*}*<sup>n</sup>* <sup>⊂</sup> *<sup>L</sup>*2(0, *<sup>T</sup>*; *<sup>L</sup>*2(Ω)).

Next, we prove part (2).

By taking *<sup>ψ</sup>* <sup>=</sup> *<sup>u</sup>*\$*<sup>m</sup> <sup>n</sup>* in (34) and by using the elementary inequality (*<sup>a</sup>* <sup>−</sup> *<sup>b</sup>*, *<sup>a</sup>*) = *<sup>a</sup>*2−*b*<sup>2</sup> <sup>2</sup> + (*a*−*b*)<sup>2</sup> <sup>2</sup> valid for all *<sup>a</sup>*, *<sup>b</sup>* <sup>∈</sup> <sup>R</sup>, we have that

$$\|\|\hat{u}\_n^m\|\|\_2^2 - \|\|\hat{u}\_n^{m-1}\|\|\_2^2 + \|\|\hat{u}\_n^m - \hat{u}\_n^{m-1}\|\|\_2^2 + \kappa\_n \nu \|\|\nabla \hat{u}\_n^m\|\|\_2^2 = 0. \tag{35}$$

Then, for any fixed *m* ∈ {1, ..., *n*} we have that

$$\|\|\hat{u}\_n^m\|\|\_2^2 \le \|\|u\_0^n\|\|\_2^2 \le \|\|u\_0\|\|\_{2^\*}^2\tag{36}$$

$$\kappa \le \sum\_{i=1}^{m} ||\nabla \hat{u}\_n^i||\_2^2 \le ||u\_0^n||\_2^2 \le ||u\_0||\_{2'}^2\tag{37}$$

$$\sum\_{i=1}^{m} ||\hat{u}\_n^i - \hat{u}\_n^{i-1}||\_2^2 \le ||u\_0^n||\_2^2 \le ||u\_0||\_2^2. \tag{38}$$

By using (36)–(38) and (33) we easily have that

{*un*}*<sup>n</sup>* is bounded in *<sup>L</sup>*∞(0, *<sup>T</sup>*; *<sup>H</sup>*), (39)

$$L^{\{\upsilon^{\eta}\}\_{\eta}}\text{ is bounded in }L^{\infty}(0,T;H)\cap L^{2}(0,T;V). \tag{40}$$

We want to prove a uniform bound in *<sup>L</sup>*2(0, *<sup>T</sup>*; *<sup>V</sup>*) also for {*un*}*n*. By a direct calculation we have that

$$\begin{split} \int\_{0}^{T} \|\nabla\boldsymbol{u}^{n}(t)\|\_{2}^{2} dt &= \sum\_{m=1}^{n} \int\_{t\_{m-1}}^{t\_{m}} \left(1 - \frac{(t - t\_{m-1})}{\kappa\_{n}}\right)^{2} \|\nabla\hat{\boldsymbol{u}}\_{n}^{m-1}\|\_{2}^{2} dt \\ &+ 2\sum\_{m=1}^{n} \int\_{t\_{m-1}}^{t\_{m}} \left(1 - \frac{(t - t\_{m-1})}{\kappa\_{n}}\right) \left(\frac{(t - t\_{m-1})}{\kappa\_{n}}\right) (\nabla\hat{\boldsymbol{u}}\_{n}^{m-1}, \nabla\hat{\boldsymbol{u}}\_{n}^{m}) dt \\ &+ \sum\_{m=1}^{n} \int\_{t\_{m-1}}^{t\_{m}} \left(\frac{(t - t\_{m-1})}{\kappa\_{n}}\right)^{2} \|\nabla\hat{\boldsymbol{u}}\_{n}^{m}\|\_{2}^{2} dt \\ &\leq \frac{\kappa\_{n}}{2} \sum\_{m=1}^{n} \|\nabla\hat{\boldsymbol{u}}\_{n}^{m-1}\|\_{2}^{2} + \frac{\kappa\_{n}}{2} \sum\_{m=1}^{n} \|\nabla\hat{\boldsymbol{u}}\_{n}^{m}\|\_{2}^{2} \\ &\leq \frac{\kappa\_{n}}{2} \|\nabla\hat{\boldsymbol{u}}\_{n}^{0}\|\_{2}^{2} + \kappa\_{n} \sum\_{m=1}^{n} \|\nabla\hat{\boldsymbol{u}}\_{n}^{m}\|\_{2}^{2}. \end{split}$$

By using (32) we obtain

$$\begin{aligned} \int\_0^T ||\nabla u^n(t)||\_2^2 \, dt &\lesssim \kappa\_n \sum\_{m=1}^n ||\nabla \hat{u}\_n^m||\_2^2 + \kappa\_n ||\nabla u\_0^n||\_2^2 \\ &\lesssim \kappa\_n \sum\_{m=1}^n ||\nabla \hat{u}\_n^m||\_2^2 + ||u\_0||\_2^2 \lesssim ||u\_0||\_2^2. \end{aligned}$$

where we have also used that *<sup>κ</sup><sup>n</sup>* <sup>=</sup> *<sup>T</sup>*/*<sup>n</sup>* and (37). Therefore we have that {*un*}*<sup>n</sup>* is bounded in *<sup>L</sup>*2(0, *<sup>T</sup>*; *<sup>V</sup>*) and then, taking into account (39), {*un*}*<sup>n</sup>* satisfies the condition *(1)* in Definition 2. Next, we show that {*un*}*<sup>n</sup>* satisfies the condition *(2)* of Definition 2. First, for all *<sup>φ</sup>* ∈ V(Ω) and *<sup>χ</sup>* <sup>∈</sup> *<sup>C</sup>*<sup>∞</sup> *<sup>c</sup>* ([0, *T*)) we have, by using (33) and (34), that

$$\int\_{0}^{T} \left( (\boldsymbol{u}^{n}(t), \boldsymbol{\phi}) \boldsymbol{\chi}(t) + ((\boldsymbol{v}^{n}(t) \cdot \nabla)\boldsymbol{v}^{n}(t), \boldsymbol{\phi}) \boldsymbol{\chi}(t) + \boldsymbol{\nu}(\nabla \boldsymbol{v}^{n}(t), \nabla \boldsymbol{\phi}) \boldsymbol{\chi}(t) \right) dt - (\boldsymbol{u}^{n}\_{0}, \boldsymbol{\phi}) \boldsymbol{\chi}(0) = 0.$$

If we define

$$R^n\_\phi := \left(\upsilon^n(t) - \mathfrak{u}^n(t), \upsilon \Delta \phi\right) + \left(\left(\upsilon^n(t) - \mathfrak{u}^n(t)\right) \otimes \upsilon^n(t) + \mathfrak{u}^n(t) \otimes \left(\upsilon^n(t) - \mathfrak{u}^n(t)\right), \nabla \phi\right),$$

then {*un*}*<sup>n</sup>* satisfies the formulation (11) and we only need to prove that (12). To this end we note that

$$\begin{split} \int\_{0}^{T} |R^{\mathrm{H}}\_{\boldsymbol{\Phi}}(t)| \, dt \leq & c \| \nabla \boldsymbol{\Phi} \|\_{\infty} \left( \int\_{0}^{T} \| \boldsymbol{u}^{\mathrm{u}}(\boldsymbol{s}) - \boldsymbol{v}^{\mathrm{u}}(\boldsymbol{s}) \|\_{2}^{2} \, ds \right)^{\frac{1}{2}} \left( \int\_{0}^{T} \| \boldsymbol{u}^{\mathrm{u}}(\boldsymbol{s}) \| ^{2} + \| \boldsymbol{v}^{\mathrm{u}}(\boldsymbol{s}) \|\_{2}^{2} \, ds \right)^{\frac{1}{2}} \\ &+ \nu \, T^{\frac{1}{2}} \| \nabla^{2} \boldsymbol{\Phi} \|\_{2} \left( \int\_{0}^{T} \| \boldsymbol{u}^{\mathrm{u}}(\boldsymbol{s}) - \boldsymbol{v}^{\mathrm{u}}(\boldsymbol{s}) \|\_{2}^{2} \, ds \right)^{\frac{1}{2}}. \end{split}$$

By a direct calculation, we have that

$$\int\_0^T ||u^n(t) - v^n(t)||\_2^2 dt = \frac{\kappa\_n}{3} \sum\_{m=1}^n ||\hat{u}\_n^m - \hat{u}\_n^{m-1}||\_2^2 \le \mathbb{C} \kappa\_n. \tag{41}$$

and therefore

$$\int\_0^T |R^n\_\Phi(t)| \, dt \le C \| |\nabla \Phi| \|\_\infty \kappa\_{n\nu} \|$$

and (12) follows.

In conclusion, we have proved that {*un*}*<sup>n</sup>* satisfies the conditions (1) and (2) of Definition 2 and thanks to Theorem 1 and Remark 3, there exists *u* satisfying the condition (1) and (2) in Definition 1. Then, in order to conclude, we only need to prove that *u* satisfies also the energy inequality. First, we note that by using (35) and (33), a direct calculation implies that for any *t* ∈ (0, *T*)

$$\|\boldsymbol{v}^{\boldsymbol{u}}(t)\|\_{2}^{2} + 2\nu \int\_{0}^{t} \|\nabla \boldsymbol{v}^{\boldsymbol{u}}(s)\|\_{2}^{2} ds \le \|\boldsymbol{u}\_{0}^{\boldsymbol{u}}\|\_{2}^{2}.\tag{42}$$

By using (40) and (41) we can infer that *v<sup>n</sup>* converges to the same limit of *un*, namely that

$$\begin{aligned} \upsilon^{\mathfrak{n}} &\rightarrow \mathfrak{u} &\text{strongly in } L^2(0, T; L^2(\Omega)),\\ \nabla \upsilon^{\mathfrak{n}} &\rightharpoonup \nabla \mathfrak{u} &\text{weakly in } L^2(0, T; L^2(\Omega)). \end{aligned} \tag{43}$$

For *<sup>k</sup>* <sup>∈</sup> <sup>N</sup> and *<sup>t</sup>* <sup>∈</sup> (0, *<sup>T</sup>*), let *<sup>χ</sup><sup>t</sup> <sup>k</sup>* be the same function already defined in (21). Noticing that <sup>−</sup>*χ*˙*<sup>t</sup> <sup>k</sup>* is positive, after multiplying (42) and integrating in time we get that, for any *t* ∈ (0, *T*), it holds

$$\frac{1}{k} \int\_{t}^{t+k} \|v^{\mathfrak{n}}(s)\|\_{2}^{2} ds + 2\nu \int\_{0}^{T} \chi\_{k}^{t}(s) \|\nabla v^{\mathfrak{n}}(s)\|\_{2}^{2} ds \leq \|u\_{0}^{\mathfrak{n}}\|\_{2}^{2} \int\_{0}^{T} (-\chi\_{k}^{t}(s)) \, ds = \|u\_{0}^{\mathfrak{n}}\|\_{2}^{2}.$$

By using (43) we get

$$\frac{1}{k} \int\_{t}^{t+k} \|\boldsymbol{u}(s)\|\_{2}^{2} \, ds + 2\nu \int\_{0}^{T} \chi\_{k}^{t}(s) \|\nabla \boldsymbol{u}(s)\|\_{2}^{2} \, ds \le \|\boldsymbol{u}\_{0}\|\_{2'}^{2}$$

and by Lebesgue differentiation and dominated convergence theorems we obtaine that for a.e. *t* ∈ (0, *T*)

$$\|\mu(t)\|\_{2}^{2} + 2\nu \int\_{0}^{t} \|\nabla \mu(s)\|\_{2}^{2} ds \le \|\mu\_{0}\|\_{2}^{2}. \tag{44}$$

Let N ⊂ (0, *T*) the set of measure zero where (44) does not hold and fix *t* ∈ N . Then, there exists {*tk*}*<sup>k</sup>* ⊂ (0, *T*) \ N such that *tk* → *t* and

$$\|\|\boldsymbol{\mu}(t\_k)\|\|\_2^2 + 2\nu \int\_0^{t\_k} \|\|\nabla \boldsymbol{\mu}(s)\|\|\_2^2 \, ds \le \|\|\boldsymbol{\mu}\_0\|\|\_2^2.$$

Since *<sup>u</sup>* <sup>∈</sup> *Cw*(0, *<sup>T</sup>*; *<sup>H</sup>*) and <sup>∇</sup>*u*(·) <sup>2</sup> <sup>2</sup> <sup>∈</sup> *<sup>L</sup>*1(0, *<sup>T</sup>*) it follows that

$$\|\|u(t)\|\|\_{2}^{2} + 2\nu \int\_{0}^{t} \|\nabla u(s)\|\|\_{2}^{2} ds \le \liminf\_{k \to \infty} \|u(t\_{k})\|\_{2}^{2} + \lim\_{k \to \infty} 2\nu \int\_{0}^{t\_{k}} \|\nabla u(s)\|\_{2}^{2} ds \le \|u\_{0}\|\_{2}^{2}.$$

and therefore (44) holds for any *t* ∈ (0, *T*).

#### *5.4. Smagorinsky-Ladyžhenskaya Model*

In this section we show how the approximation by adding a nonlinear stress tensor produce weak solutions. We consider for *<sup>n</sup>* <sup>∈</sup> <sup>N</sup> the following boundary initial value problem

$$\begin{cases} \partial\_t \, u^n + \left( u^n \cdot \nabla \right) u^n + \nabla p^n - \nu \Delta u^n - \frac{1}{n} \operatorname{div} (|Du^n| D u^n) = 0 & \text{in } (0, T) \times \Omega, \\\\ \operatorname{div} \, u^n = 0 & \text{in } (0, T) \times \Omega, \\ u^n = 0 & \text{on } (0, T) \times \partial \Omega, \\ u^n(0) = u\_0 & \text{in } \Omega, \end{cases} \tag{45}$$

where *Du<sup>n</sup>* <sup>=</sup> <sup>∇</sup>*un*+(∇*un*)*<sup>T</sup>* <sup>2</sup> . This system has been introduced for numerical approximation of turbulent flows by Smagorinsky [22] and its analysis as a possible approximation for the Navier-Stokes equations started with the studies by Ladyženskaya [23], cf. also Reference [24] for the role of this method in the analysis of Large Eddy Simulation models. For the analysis also of related models, with general stress tensor given by *S*(*v*) = *S*(*Dv*) = |*Dv*| *<sup>p</sup>*−2*Dv*, with various values of *p*, see References [18,25] and also the more recent Reference [26].

**Theorem 5.** *Let u*<sup>0</sup> ∈ *H. Then, it holds that*


$$
\mu^n \to \mu \qquad \text{strongly in } L^2(0, T; L^2(\Omega)).
$$

**Proof.** By using the theory of monotone operators (cf. References [8,18]) there exists a unique *<sup>u</sup><sup>n</sup>* <sup>∈</sup> *<sup>C</sup>*(0, *<sup>T</sup>*; *<sup>L</sup>*2(Ω)) weak solution of (45) with *<sup>u</sup><sup>n</sup>* <sup>∈</sup> *<sup>L</sup>*∞(0, *<sup>T</sup>*; *<sup>H</sup>*) <sup>∩</sup> *<sup>L</sup>*2(0, *<sup>T</sup>*; *<sup>V</sup>*) and *Du<sup>n</sup>* <sup>∈</sup> *<sup>L</sup>*3(0, *<sup>T</sup>*; *<sup>L</sup>*3(Ω)) such that it holds

$$\|\|u^n(t)\|\|\_{2}^{2} + 2\nu \int\_{0}^{t} \|\nabla u^n(s)\|\|\_{2}^{2} \, ds + \frac{2}{n} \int\_{0}^{t} \|\|Du^n(s)\|\|\_{3}^{3} \, ds \le \|u\_0^n\|\|\_{2}^{2}.$$

Observe that by Korn inequality *Du<sup>n</sup>* <sup>3</sup> ∼ 
∇*u<sup>n</sup>* 3.

To prove the second part of Theorem <sup>5</sup> we show that {*un*}*<sup>n</sup>* satisfy the conditions in Definition 2. Define the remainder *R<sup>n</sup> <sup>φ</sup>*(*t*) by

$$R^{\eta}\_{\Phi}(t) := -\frac{1}{n} \int\_{\Omega} |Du^{\eta}(t)| Du^{\eta}(t) \cdot D\phi \,d\mathbf{x}.$$

By means of the Hölder inequality we get

$$|R^n\_{\phi}(t)| \le \frac{1}{n} \int\_{\Omega} |Du^n(t)|^2 |D\phi| \, d\mathfrak{x} \le \frac{1}{n} ||Du^n||\_3^2 ||D\phi||\_3.$$

Consequently, it also holds

$$\int\_0^T |R^n\_\phi(t)| \, dt \le \frac{T^{1/3}}{n^{1/3}} \left( \frac{1}{n} \int\_0^T ||Du^n||\_3^3 \, dt \right)^{2/3} ||D\phi||\_{3\nu}.$$

showing that *R<sup>n</sup> <sup>φ</sup>* <sup>→</sup> 0 in *<sup>L</sup>*1(0, *<sup>T</sup>*). Since the other conditions in Definition <sup>2</sup> are trivially satisfied, an application of Theorem 1 finally ends the proof.

**Author Contributions:** Both authors contributed equally to Conceptualization, Writing—review & editing. Both authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** LCB and SS are member of the group GNAMPA of INdAM.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Fluids* Editorial Office E-mail: fluids@mdpi.com www.mdpi.com/journal/fluids

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-1999-9