Long-Time Reynolds Averaging of Reduced Order Models for Fluid Flows: Preliminary Results

TL;DR

This paper investigates the long-term average energy transfer among modes in reduced order models of fluid flows, combining theoretical proofs and numerical experiments using the Burgers equation.

Contribution

It provides the first analytical results on energy exchange in ROM modes and compares eigenfunction and POD modes, supported by numerical evidence.

Findings

Time-average energy transfer from low to high POD modes is positive over long intervals.

Analytical differences between eigenfunction and POD modes are highlighted.

Numerical results support theoretical predictions about energy transfer.

Abstract

We perform a theoretical and numerical investigation of the time-average of energy exchange among modes of Reduced Order Models (ROMs) of fluid flows. We are interested in the statistical equilibrium problem, and especially in the possible forward and backward average transfer of energy among ROM basis functions (modes). We consider two types of ROM modes: eigenfunctions of the Stokes operator and Proper Orthogonal Decomposition (POD) modes. We prove analytical results for both types of ROM modes and we highlight the differences between them. We also investigate numerically whether the time-average energy exchange between POD modes is positive. To this end, we utilize the one-dimensional Burgers equation as a simplified mathematical model, which is commonly used in ROM tests. The main conclusion of our numerical study is that, for long enough time intervals, the time-average energy exchange from low index POD modes to high index POD modes is positive, as predicted by our theoretical results.

Figures

Tables

Table 1. Table 1. Case 1: Time-averages of e m ​ ( u ) subscript e 𝑚 𝑢 \mathrm{e}_{m}(u) and ℰ m ​ ( u ) subscript ℰ 𝑚 𝑢 \mathcal{E}_{m}(u) for d = 37 𝑑 37 d=37 and different m 𝑚 m values.

$m$	${\int }_0^1{e}_m\left(u\left(s\right)\right)𝑑s$	${\int }_0^1{ℰ}_m\left(u\left(s\right)\right)𝑑s$
3	2.6170e-02	1.0451e-03
6	7.3208e-03	2.2524e-03
9	1.4934e-03	1.5977e-03
15	3.6181e-05	3.6287e-05
20	1.2776e-06	7.0356e-07
25	4.6207e-08	9.5638e-09
30	2.0257e-09	-1.1127e-10
35	8.2806e-11	-2.2595e-11

Table 2. Table 2. Case 2: Time-averages of e m ​ ( u ) subscript e 𝑚 𝑢 \mathrm{e}_{m}(u) and ℰ m ​ ( u ) subscript ℰ 𝑚 𝑢 \mathcal{E}_{m}(u) for d = 38 𝑑 38 d=38 , Δ ​ t = 10 − 2 Δ 𝑡 superscript 10 2 \Delta t=10^{-2} , 1000 1000 1000 equally spaced quadrature points, and different m 𝑚 m values.

$m$	$\frac{1}{10}{\int }_0^{10}{e}_m\left(u\left(s\right)\right)𝑑s$	$\frac{1}{10}{\int }_0^{10}{ℰ}_m\left(u\left(s\right)\right)𝑑s$
3	3.3652e-03	8.6523e-05
6	1.0001e-03	2.1570e-04
9	2.1132e-04	1.8614e-04
15	5.4224e-06	5.5724e-06
20	4.3119e-07	2.5667e-07
25	5.6884e-08	1.0451e-08
30	1.1327e-09	3.2736e-10
35	2.0469e-11	2.6396e-12

Table 3. Table 4. Case 2: Time-averages of e m ​ ( u ) subscript e 𝑚 𝑢 \mathrm{e}_{m}(u) and ℰ m ​ ( u ) subscript ℰ 𝑚 𝑢 \mathcal{E}_{m}(u) for d = 43 𝑑 43 d=43 , Δ ​ t = 10 − 4 Δ 𝑡 superscript 10 4 \Delta t=10^{-4} , 10000 10000 10000 equally spaced quadrature points, and different m 𝑚 m values.

$m$	$\frac{1}{10}{\int }_0^{10}{e}_m\left(u\left(s\right)\right)𝑑s$	$\frac{1}{10}{\int }_0^{10}{ℰ}_m\left(u\left(s\right)\right)𝑑s$
3	3.5637e-03	2.7841e-06
6	1.1664e-03	8.1109e-06
9	2.7346e-04	1.3956e-05
15	1.2084e-05	1.1937e-05
20	1.1785e-06	1.6370e-06
25	2.2727e-07	1.2913e-07
30	6.2606e-08	9.5191e-09
35	2.1106e-09	6.1538e-10
40	1.0330e-10	2.0241e-11

Table 4. Table 6. Case 3: Time-averages of e m ​ ( u ) subscript e 𝑚 𝑢 \mathrm{e}_{m}(u) for d = 36 𝑑 36 d=36 , Δ ​ t = 10 − 2 Δ 𝑡 superscript 10 2 \Delta t=10^{-2} , different m 𝑚 m values, and all subintervals used in the composite trapezoidal rule.

$m$	$\frac{1}{100}{\int }_0^{100}{e}_m\left(u\left(s\right)\right)𝑑s$	$\frac{1}{10}{\int }_0^{10}{e}_m\left(u\left(s\right)\right)𝑑s$	${\int }_0^1{e}_m\left(u\left(s\right)\right)𝑑s$
3	3.3687e-04	3.3683e-04	1.7634e-04
6	1.0001e-04	1.0001e-04	7.3142e-05
9	2.1146e-05	2.1147e-05	1.6701e-05
15	5.5621e-07	5.5742e-07	4.7129e-07
20	7.9605e-08	7.9605e-08	6.3799e-08
25	8.3483e-09	8.3489e-09	5.7426e-09
30	1.0839e-10	1.0840e-10	9.2430e-11
35	2.7710e-12	2.7718e-12	2.9575e-12

Table 5. Table 7. Case 3: Time-averages of ℰ m ​ ( u ) subscript ℰ 𝑚 𝑢 \mathcal{E}_{m}(u) for d = 36 𝑑 36 d=36 , Δ ​ t = 10 − 2 Δ 𝑡 superscript 10 2 \Delta t=10^{-2} , different m 𝑚 m values, and all subintervals used in the composite trapezoidal rule.

$m$	$\frac{1}{100}{\int }_0^{100}{ℰ}_m\left(u\left(s\right)\right)𝑑s$	$\frac{1}{10}{\int }_0^{10}{ℰ}_m\left(u\left(s\right)\right)𝑑s$	${\int }_0^1{ℰ}_m\left(u\left(s\right)\right)𝑑s$
3	8.6270e-06	7.0320e-06	8.1057e-05
6	2.1568e-05	2.1520e-05	3.2405e-05
9	1.8614e-05	1.8599e-05	2.0255e-05
15	5.5718e-07	5.4232e-07	5.8207e-07
20	5.6044e-08	5.5874e-08	6.1423e-08
25	1.0936e-09	9.9763e-10	5.9386e-10
30	3.3046e-11	3.2478e-11	4.0949e-11
35	6.4170e-13	6.3739e-13	1.2125e-12

Table 6. Table 8. Case 4: Time-averages of e m ​ ( u ) subscript e 𝑚 𝑢 \mathrm{e}_{m}(u) and ℰ m ​ ( u ) subscript ℰ 𝑚 𝑢 \mathcal{E}_{m}(u) for Δ ​ t = 2 ∗ 10 − 5 Δ 𝑡 2 superscript 10 5 \Delta t=2*10^{-5} , different m 𝑚 m values, d = 18 𝑑 18 d=18 , and 5000 5000 5000 equally spaced subintervals used in the composite trapezoidal rule.

$m$	$\frac{1}{0.1}{\int }_0^{0.1}{e}_m\left(u\left(s\right)\right)𝑑s$	$\frac{1}{0.1}{\int }_0^{0.1}{ℰ}_m\left(u\left(s\right)\right)𝑑s$
3	-6.8687e-04	2.7378e-05
5	-2.6333e-05	1.7795e-05
7	-9.1458e-07	4.2432e-06
9	-9.8188e-09	4.1220e-07
13	2.3800e-10	1.5749e-09
15	8.5423e-12	4.8043e-11

Equations

$$\{\begin{array}{rl}{\partial }_t\mathbf{u}-\nu \Delta \mathbf{u}+(\mathbf{u}\cdot \nabla )\mathbf{u}+\nabla p& =\mathbf{f}\text{mbox}in\Omega \times (0,T),\\ \nabla \cdot \mathbf{u}& =0\text{mbox}in\Omega \times (0,T),\\ \mathbf{u}& =0\text{mbox}on\Gamma \times (0,T),\\ \mathbf{u}(\cdot ,0)& =\mathbf{u}(0)\text{mbox}in\Omega ,\end{array}$$

$$Re=\frac{UL}{\nu },$$

$$\mathbf{u}=\sum _{k=1}^{\infty }{\mathbf{u}}_k{\mathbf{w}}_k,$$

$${\mathbf{u}}_k=\frac{{\int }_{\Omega }\mathbf{u}\cdot {\mathbf{w}}_{k}dx}{{\int }_{\Omega }|{\mathbf{w}}_k{|}^{2}dx}.$$

$$\frac{1}{2}\Vert \mathbf{u}(t){\Vert }^2+\nu {\int }_0^t\Vert \nabla \mathbf{u}(s){\Vert }^{2}ds=\frac{1}{2}\Vert \mathbf{u}(0){\Vert }^2+{\int }_0^t{\int }_{\Omega }\mathbf{f}\cdot \mathbf{u}dxds.$$

$${\partial }_t\frac{|\mathbf{u}{|}^2}{2}+|\nabla \mathbf{u}{|}^2-\Delta \frac{|\mathbf{u}{|}^2}{2}+\text{div }\left(\frac{\mathbf{u}|\mathbf{u}{|}^2}{2}+p\mathbf{u}\right)=\mathbf{f}\cdot \mathbf{u}.$$

$${\int }_0^T{\int }_{\Omega }|\nabla \mathbf{u}{|}^2\varphi dxdt={\int }_0^T{\int }_{\Omega }\left[\frac{|\mathbf{u}{|}^2}{2}\left({\partial }_t\varphi +\Delta \varphi \right)+\left(\frac{|\mathbf{u}{|}^2}{2}+p\right)\mathbf{u}\cdot \nabla \varphi +\mathbf{f}\cdot \mathbf{u}\varphi \right]dxdt$$

$$\mathrm{E}(\mathbf{u})=\frac{1}{2}\sum _{k=1}^{\infty }\Vert {\mathbf{u}}_k{\Vert }^2.$$

$${\mathbf{u}}_{m^{\prime },m^{\prime \prime }}:=\sum _{k=m^{\prime }}^{m^{\prime \prime }}{\mathbf{u}}_k{\mathbf{w}}_k,$$

$$\mathrm{E}({\mathbf{u}}_{m^{\prime },m^{\prime \prime }})=\frac{1}{2}\sum _{k=m^{\prime }}^{m^{\prime \prime }}\Vert {\mathbf{u}}_k{\Vert }^2.$$

$$\dot{a}=Aa+a^{\top }Ba,$$

$$\mathbf{u}=\mathbf{y}+\mathbf{z},$$

$$\mathbf{y}=\sum _{k=1}^m{\mathbf{u}}_k{\mathbf{w}}_k={\mathbf{P}}_m\mathbf{u}\text{and}\mathbf{z}=\sum _{k=m+1}^{+\infty }{\mathbf{u}}_k{\mathbf{w}}_k=(\mathbf{I}-{\mathbf{P}}_m)\mathbf{u}=:{\mathbf{Q}}_m\mathbf{u}.$$

$$\mathcal{V}=\{\mathbf{\phi }\in \mathcal{D}(\Omega {)}^3,\nabla \cdot \mathbf{\phi }=0\},$$

$$H=\left\{\mathbf{u}\in L^2(\Omega {)}^3,\nabla \cdot \mathbf{u}=0,\mathbf{u}\cdot \mathbf{n}=0\text{ on }\Gamma \right\},$$

$$V=\left\{\mathbf{u}\in H_0^1(\Omega {)}^3,\nabla \cdot \mathbf{u}=0\right\},$$

$$\frac{d\mathbf{u}}{dt}+\nu A\mathbf{u}+B(\mathbf{u},\mathbf{u})=\mathrm{P}\mathbf{f},$$

$${\mathbf{P}}_m\mathbf{f}=\mathbf{f}.$$

$$\frac{d\mathbf{y}}{dt}-\nu {\mathbf{P}}_m(\Delta \mathbf{u})+{\mathbf{P}}_{m}B(\mathbf{y}+\mathbf{z},\mathbf{y}+\mathbf{z})={\mathbf{P}}_m\mathbf{f},$$

$$\frac{d\mathbf{z}}{dt}-\nu {\mathbf{Q}}_m(\Delta \mathbf{u})+{\mathbf{Q}}_{m}B(\mathbf{y}+\mathbf{z},\mathbf{y}+\mathbf{z})=0,$$

$$-\nu {\int }_{\Omega }{\mathbf{P}}_m(\Delta \mathbf{u})\cdot \mathbf{y}dx=-\nu {\int }_{\Omega }(\Delta \mathbf{u})\cdot \mathbf{y}dx=-\nu {\int }_{\Omega }(\Delta \mathbf{y}+\Delta \mathbf{z})\cdot \mathbf{y}dx=\nu \Vert \nabla \mathbf{y}{\Vert }^2+\nu {\int }_{\Omega }\nabla \mathbf{y}:\nabla \mathbf{z}dx,$$

$$-\nu {\int }_{\Omega }{\mathbf{Q}}_m(\Delta \mathbf{u})\cdot \mathbf{z}dx=-\nu {\int }_{\Omega }(\Delta \mathbf{u})\cdot \mathbf{z}dx=-\nu {\int }_{\Omega }(\Delta \mathbf{y}+\Delta \mathbf{z})\cdot \mathbf{z}dx=\nu \Vert \nabla \mathbf{z}{\Vert }^2+\nu {\int }_{\Omega }\nabla \mathbf{y}:\nabla \mathbf{z}dx.$$

$$\frac{1}{2}\frac{d}{dt}\Vert \mathbf{y}{\Vert }^2+\nu \Vert \nabla \mathbf{y}{\Vert }^2+\nu (\nabla \mathbf{y},\nabla \mathbf{z})+(B(\mathbf{y}+\mathbf{z},\mathbf{y}+\mathbf{z}),\mathbf{y})=(\mathbf{f},\mathbf{y}),$$

$$\frac{1}{2}\frac{d}{dt}\Vert \mathbf{z}{\Vert }^2+\nu \Vert \nabla \mathbf{z}{\Vert }^2+\nu (\nabla \mathbf{y},\nabla \mathbf{z})+(B(\mathbf{y}+\mathbf{z},\mathbf{y}+\mathbf{z}),\mathbf{z})\le 0,$$

$$\frac{1}{2}\frac{d}{dt}\mathrm{E}(\mathbf{y})+\nu \Vert \nabla \mathbf{y}{\Vert }^2+\nu (\nabla \mathbf{y},\nabla \mathbf{z})=(B(\mathbf{y},\mathbf{y}),\mathbf{z})-(B(\mathbf{z},\mathbf{z}),\mathbf{y})+(\mathbf{f},\mathbf{y}),$$

$$\frac{1}{2}\frac{d}{dt}\mathrm{E}(\mathbf{z})+\nu \Vert \nabla \mathbf{z}{\Vert }^2+\nu (\nabla \mathbf{y},\nabla \mathbf{z})\le -(B(\mathbf{y},\mathbf{y}),\mathbf{z})+(B(\mathbf{z},\mathbf{z}),\mathbf{y}).$$

$${\mathrm{e}}_m(\mathbf{u}):=e^{\uparrow }(\mathbf{u})-e^{\downarrow }(\mathbf{u})\text{ with}{e}^{\uparrow }(\mathbf{u}):=-(B(\mathbf{y},\mathbf{y}),\mathbf{z}),e^{\downarrow }(\mathbf{u}):=(B(\mathbf{z},\mathbf{z}),\mathbf{y}).$$

$${\mathcal{E}}_m(\mathbf{u}):=-\nu (\nabla \mathbf{y},\nabla \mathbf{z}).$$

$$\frac{1}{2}\frac{d}{dt}\mathrm{E}(\mathbf{y})+\nu \Vert \nabla \mathbf{y}{\Vert }^2={\mathcal{E}}_m(\mathbf{u})-{\mathrm{e}}_m(\mathbf{u})+(\mathbf{f},\mathbf{y}),$$

$$\frac{1}{2}\frac{d}{dt}\mathrm{E}(\mathbf{z})+\nu \Vert \nabla \mathbf{z}{\Vert }^2\le {\mathcal{E}}_m(\mathbf{u})+{\mathrm{e}}_m(\mathbf{u}).$$

$$\nu {\int }_{\Omega }\nabla \mathbf{y}:\nabla \mathbf{z}dx,$$

$$A{\mathcal{W}}_k={\lambda }_k{\mathcal{W}}_k\text{and}{\int }_{\Omega }{\mathcal{W}}_k\cdot {\mathcal{W}}_{j}dx={\delta }_{kj}.$$

$${\int }_{\Omega }\nabla {\mathcal{W}}_k:\nabla {\mathcal{W}}_{j}dx=0\text{for }k\ne j,$$

$$\mathbf{u}=\mathbf{y}+\mathbf{z}:=\sum _{k=1}^m{c}_k{\mathcal{W}}_k+\sum _{k=m+1}^{\infty }{c}_k{\mathcal{W}}_k={\mathbf{P}}_m\mathbf{u}+{\mathbf{Q}}_m\mathbf{u},$$

Full text

1. Introduction

In this note we combine some results on the long-time averaging of fluid equations with the recently developed techniques for reduced order model (ROM) development. In this preliminary work we start proving some analytical results that characterize the time-averaged effect of the exchange of energy between various modes, both in the case of the computable decomposition made with proper orthogonal decomposition (POD) type basis functions and with the abstract basis made with eigenfunctions. We will show that the results obtainable with a generic (but computable) basis are less precise than those obtainable with the abstract spectral basis, the difference coming from the lack of orthogonality of the gradients of the POD basis functions.

We then provide a few numerical examples. Concerning the analytical results we will prove partial results for the energy exchange between large and small scales, showing the difference between the use of a spectral type basis, versus a POD one. In particular, we are interested in results connected to the statistical equilibrium problem, which can be deduced in a computable way by a long-time averaging of the solutions. We want to investigate the possible forward and backward average transfer of energy. The properties of a turbulent flow are computable (and relevant) only in an average sense. In this respect, we want to follow the most classical approach dating back to Stokes, Reynolds, and Prandtl of considering long-time averages of the solution as the key quantity to be computed or observed. Therefore, we do not need to consider statistical averages and link them with time averaging by means of (unproved) ergodic hypotheses.

To introduce the problem that we will consider, we recall that a Newtonian incompressible flow (with constant density) can be described by the Navier-Stokes equations (NSE in the sequel)

[TABLE]

with Dirichlet conditions when the motion takes place in a smooth and bounded domain $\Omega\subset\mathbbm{R}^{3}$ with solid walls $\Gamma:=\partial\Omega$. The unknowns are the velocity field ${\bf u}$ and the scalar pressure $p$, while the positive constant $\nu>0$ is the kinematic viscosity. The key parameter to detect the nature of the problem is the non-dimensional Reynolds number, which is defined as

[TABLE]

where $U$ and $L$ are a characteristic velocity and length of the problem. In realistic problems, the Reynolds number can be extremely large (in many cases of the order of $10^{8}$, but up to the order of $10^{20}$ in certain geophysical problems). For simplicity in the notation, we use as a control parameter the viscosity and we assume that it is very small. Hence, the effect of the regularization (similar to the diffusion in heat transfer) due to the Laplacian is negligible and the behavior of solutions is really turbulent and rather close to the motion of ideal fluids. Due to the well-known difficulties in performing direct numerical simulations (DNS), it is nowadays a well-established technique that of trying to reduce the computational efforts by simulating only the largest scales, which are nevertheless the only ones really observable and the only ones which are needed in order to optimize macroscopic properties as, e.g., drag or lift in the craft design. In this framework, the large eddy simulation (LES) methods, which emerged in the last 30 years, are among the most popular, and they found a very relevant role within both theorists’ and practitioners’ communities. For recent LES reviews see for instance the monographs [31, 2, 20, 5].

The LES methods are in many cases very well developed and both theoretically and computationally appealing, especially for problems without boundaries, but many difficulties and open problems arise when facing solid boundaries. In most cases the design of efficient LES methods is guided by deep properties of the solutions, as emerging from fine theorems of mathematical analysis. Furthermore, the ultimate goal of having a golden standard is far from being obtained, and large families of methods (wave-number asymptotics, differential filters, $\alpha$-models) attracted the interest of different communities, spanning from the pure mathematicians, to the applied geophysicist and mechanical engineers, as well as computational practitioners. For this reasons we believe that it is important to have some well-defined and clearly stated guidelines, that can be adapted to different problems. In this way the methods can be improved with insight not only from experts in modeling, but also from mathematicians, physicists, and computational scientists.

In this respect, we point out that very recently the use of other (more flexible and computationally simpler) ways of finding approximate systems has become very popular. The LES methods itself can be specialized or even glued with other ways of determining approximate and much smaller systems, which are computable in a very efficient way. For instance, reduced order modeling is increasingly becoming an accepted paradigm, in which applications to fluids are still being developed [14, 18, 27, 28, 30].

The basic ansatz at the basis of the use of these models is the approximation of the velocity by a truncation of the series

[TABLE]

where $\{{\bf w}_{k}\}_{k\in\mathbbm{N}}$ is a basis constructed by using the POD, not necessarily made with eigenfunctions of the Stokes operator, and the coefficients of the $L^{2}$-projections are evaluated as follows

[TABLE]

The appealing property of this approach is that the choice of the basis is adapted-to and determined-by the problem itself. Generally the basis is chosen after a preliminary numerical computation, hence it contains at least the basic features of the solution and geometry of the problem to be studied. The other basic fact is that the kinetic energy of the problem is the key quantity under consideration; in fact the $L^{2}$-projection is used to determine the approximate velocity and also the energy content in the basis construction. To determine the number $r\in\mathbbm{N}$ such that the solution is projected over the span generated by the (orthogonal) functions $V_{r}:=\text{span}\{{\bf w}_{1},\dots,{\bf w}_{r}\}$, generally it is assumed that the projection of the solution over $V_{r}$ contains a large percentage (say 80$\%$) of the total kinetic energy of the underlying system.

It turns out that a basis associated with the problem at hand can greatly improve the effectiveness of the ROM. Its proper choice can be of great interest in applications to fluid flows [32, 33, 35]. The main goal and novelty of the present paper is to try to combine results on the long-time behavior of fluid flows (especially in the case of statistical equilibrium) with reduced order modeling. We use this approach in order to capture both the long-time averages and the energy exchange of the ROM solutions. In this respect, we are extending to the POD setting the results based on statistical solutions by Foias et al. [10, 11] and those more recently obtained for time-averages in [3, 21]. In this respect, the main theoretical results of this paper, stated in Theorems 4.2 and 4.3 below, can be viewed as a spectral version of the results of [3, 21]. It tells that low frequency modes yield a Reynolds stress which is dissipative in the mean, the total spatial mean work of it being larger than the long time average of the dissipation of the fluctuations, which is consistent with observations and results in [3, 10]. However, the analysis shows that the triad interaction between high and low frequency modes yields an additional non positive term in the budget between the Reynolds stress of high modes and the corresponding mean dissipation. This term may be non dissipative and may permit an inverse energy cascade, which is not in contradiction with the fact that the total Reynolds stress is dissipative in mean.

Plan of the paper. The paper is organized as follows: In Section 2, we outline the general framework for ROMs of fluid flows, and we display the exchange of energy between large scales and small scales for two ROM bases: the POD and the Stokes eigenfunctions. In Section 3, we present some preliminaries on long-time averages. In Section 4, we prove the main theoretical results for the average transfer of energy for ROMs constructed with the POD and the Stokes eigenfunctions. In Section 5, we investigate the theoretical results in the numerical simulation of the one-dimensional Burgers equation. Finally, in Section 6, we draw conclusions and outline future research directions.

2. Reduced Order Modeling

As outlined in the introduction, one key quantity in the pure and applied analysis of the Navier-Stokes equations is the kinetic energy, since it is both a meaningful physical quantity, but also the analysis of its budget is at the basis of the abstract existence results (cf. Constantin and Foias [6]) and also of the conventional turbulence theories of Kolmogorov [17]. It is well-known that after testing the NSE (1.1) by ${\bf u}$ and integrating over the space-time, one (formally) obtains the global energy balance

[TABLE]

At present, we are only able to prove that the above balance holds true with the sign of “less or equal” for the class of weak (or turbulent) solutions that we are able to construct globally in time, without restrictions on the viscosity and on the size of square summable initial velocity ${\bf u}(0)$ and external force ${\bf f}$. It is of fundamental importance in many problems in pure mathematics to understand under which hypotheses the equality holds true. We are now focusing on the “global energy” which is an averaged quantity, since it is the integral of the square modulus of the velocity over the entire domain. We also point out that at the other extreme one can deduce, without the integration over the domain, the point-wise relation

[TABLE]

In between there is the so called “local energy” which can be obtained by multiplying the NSE by ${\bf u}\,\phi$, where $\phi$ is a bump function, before integrating in the space-time variables. The goal is to show that

[TABLE]

holds (at least with the inequality sign) for all smooth scalar functions $\phi\in C^{\infty}_{0}((0,T)\times\Omega)$ such that $\phi\geq 0$. The validity of such an inequality is one of the requests to use the partial regularity results, but it is also one of the requests to be satisfied by weak solutions constructed by numerical or LES methods. In this respect, see Guermond, Oden and Prudhomme [13] and also [4].

In this paper we study the global energy in the perspective that it can be reconstructed in a computable way or it can be well approximated by the POD basis functions $\{{\bf w}_{k}\}$.

The fact that the functions $\{{\bf w}_{k}\}$ can be constructed to be orthonormal with respect to the scalar product in $(L^{2}(\Omega)^{3},\,\|\cdot\|)$ allows us to evaluate the kinetic energy easily by the following numerical series

[TABLE]

Since we are going to use only a reduced number of ROM modes, it is relevant to consider the energy contained in functions described by a restricted set of indices. Hence, following the notation in [11], if we define

[TABLE]

then the kinetic energy content of ${\bf u}_{m^{\prime},m^{\prime\prime}}$ is simply evaluated as follows

[TABLE]

We denote by $\mathbf{P}_{m}$ the projection operator over the subspace $V_{m}$ spanned by the first $m$-functions $\{{\bf w}_{k}\}_{k=1,\dots,m}$ and we want to investigate the energy transfer between the various modes, together with averaged long-time behavior associated with this splitting.

We are going to adapt well-known studies on the decomposition in small and large eddies. This would be the case if the functions ${\bf w}_{k}$ are chosen to be the eigenfunctions of the Stokes operator, hence associated with large and small frequencies. In our setting the basis is determined by the solution of a simplified problem, which can be treated computationally, and the basis functions are orthonormal in $L^{2}(\Omega)$, but we cannot expect that their gradients are also orthogonal.

For the NSE, the standard ROM is constructed as follows:

1. (i)

choose modes $\{{\bf w}_{1},\ldots,{\bf w}_{d}\}$, which represent the recurrent spatial structures of the given flow; 2. (ii)

choose the dominant modes $\{{\bf w}_{1},\ldots,{\bf w}_{m}\}$, with $m\leq d$, as basis functions for the ROM; 3. (iii)

use a Galerkin truncation ${\bf u}_{m}=\sum_{j=1}^{m}a_{j}\,{\bf w}_{j}$; 4. (iv)

replace ${\bf u}$ with ${\bf u}_{m}$ in the NSE; 5. (v)

use a Galerkin projection of NSE (${\bf u}_{m}$) onto the ROM space $V_{m}:=\mbox{span}\{{\bf w}_{1},\ldots,{\bf w}_{m}\}$ to obtain a low-dimensional dynamical system, which represents the ROM:

[TABLE]

where $a$ is the vector of unknown ROM coefficients and $A,B$ are ROM operators; 6. (vi)

in an offline stage, compute the ROM operators; 7. (vii)

in an on-line stage, repeatedly use the ROM (for various parameter settings and/or longer time intervals).

Hence, there is a very natural splitting of the velocity field ${\bf u}$ into two components, the part coherent with the basis expansion associated with the less energetic modes and the remainder. This can be formalized as follows:

[TABLE]

where, for $m\in\mathbbm{N}$ which can be selected computationally (based, e.g., on the relative kinetic energy content in the first $m$ POD-modes, but other choices relative to the enstrophy are possible) in order to have a significant amount of the energy content of the flow

[TABLE]

We observe that we are considering the functions $\{{\bf w}_{k}\}$ as divergence-free. Generally they are not “exactly divergence-free,” but numerically we can consider that they have vanishing divergence, hence in the computations which will follow the pressure terms can be dropped by a standard Leray projection. It will be nevertheless interesting to consider also bases which are not divergence-free. Relevant results, due to the computational simplifications, could be also derived by the combination of ROM with artificial compressibility methods, as those introduced in [8, 12].

In addition, we consider the external force as stationary, that is ${\bf f}={\bf f}(x)\in L^{2}(\Omega)^{3}$ and we look for conditions holding at statistical equilibrium. Our purpose is to determine –if possible– the long-time behavior of ${\bf y}$ and to analyze the energy budget between low and high modes in the orthogonal decomposition determined by the functions ${\bf w}_{k}$.

As usual in many problems fluid mechanics, we use the Hilbert space functional setting with

[TABLE]

where $\mathbf{n}$ denotes the outward normal unit vector. Moreover, $V^{\prime}$ is the topological dual space of $V$. We will also denote by $<\,.\,,\,.\,>$ the duality pairing between $V$ and $V^{\prime}$. We recall that $\mathcal{V}$ is dense in $H$ and $V$ for their respective topologies [9, 22].

Once we project $L^{2}(\Omega)^{3}$ over the subspace $H$ of divergence-free and tangential vector fields by means of the Leray projection operator $\mathrm{P}$, we have the following abstract (functional) equation in $H$

[TABLE]

where $A:=-\mathrm{P}\,\Delta$, while $B({\bf u},{\bf u}):=\mathrm{P}\,(({\bf u}\cdot\nabla)\,{\bf u})$. As usual in this analysis (see for instance [11]) we can start by assuming that the input force can be decomposed within a finite sum of basis functions (or that it belongs to $V_{m}$, which will be clarified in the next section section, in particular by Theorem 2.1), hence

[TABLE]

We split the Navier-Stokes equations into a coupled system for ${\bf y}\in\mathbf{P}_{m}H$ and ${\bf z}\in(\mathbf{P}_{m}H)^{\perp}=\mathbf{Q}_{m}H$ as follows

[TABLE]

where we used that both $\mathbf{P}_{m}$ and $\mathbf{Q}_{m}$ commute with the time derivative.

Once we evaluate the kinetic energy, since $\mathbf{P}_{m}{\bf y}={\bf y}$ and $\mathbf{Q}_{m}\,{\bf z}={\bf z}$ we get (by integrating by parts and by using the fact that functions vanish at the boundary) that

[TABLE]

In this way we can obtain the following equality and inequality

[TABLE]

These are the two basic balance equations that we will use to infer the behavior and transfer of the kinetic energy between ${\bf y}$ and ${\bf z}$. Notice that the balance relation for ${\bf y}$, involving just a finite combination of rather smooth functions is an equality, while the second one is an inequality. In fact, the second one can be derived by a limiting argument and in the limit the lower semi-continuity of the norm will produce the inequality.

Since the tri-linear term $(B({\bf u},{\bf u}),{\bf w})$ is skew-symmetric with respect to the last two variables, we obtain from (2.2)

[TABLE]

This is a formal setting, which is obviously true for strong solutions of the NSE, where the inequality in (2.3) is an equality. When considering weak solutions, the integral $(B({\bf z},{\bf z}),{\bf z})$ might be not defined in $L^{1}(0,T)$ for regularity issues. However, one can still rigorously derive (2.3) by a double frequency truncation or a regularization of the operator $B$ by considering $(B({\bf z}\star\rho_{\varepsilon},{\bf z}),{\bf z})$ for a given standard mollifier $\rho_{\varepsilon}$ and passing to the limit when $\varepsilon\to 0$. Details are standard and out of the scope of the present paper.

We observe that $-(B({\bf y},{\bf y}),{\bf z})$ is the energy flux induced in the more energetic terms by the inertial forces associated to less energetic modes, while $-(B({\bf z},{\bf z}),{\bf y})$ is the energy flux induced in the less energetic terms by the inertial forces associated to more energetic modes. In a schematic way we can decompose the rate of transfer of kinetic energy $\mathrm{e}_{m}({\bf u})$ into two terms as follows

[TABLE]

We also use the following notation:

[TABLE]

Hence, we can rewrite (2.3) as follows

[TABLE]

Remark 2.1**.**

We recall that apart from extremely simple geometries and provided one is willing to use in a systematic way special functions as the Bessel ones or the spherical harmonics (which are nevertheless time consuming in their evaluation), the explicit calculations in numerical tests will not be so easy to be obtained in a precise and efficient way. Hence, the solution of (2.1) and the long-time integration of its solutions poses hard numerical problems.

We observe that, since the viscosity $\nu>0$ can be considered small in applications to real-life flows, and since the interaction between the more energetic and less energetic basis functions is weak, we can suppose that the integral

[TABLE]

even if non vanishing is negligible in the long-time (at least in a first approximation). We will show later on that –at least on the numerical tests– this assumption is reasonable and that the results obtained by using this approximation are numerically sound.

We point out for the reader that it we have a first fundamental difference with respect to the classical splitting based on the use of a spectral basis (which will be recalled in Section 2.1, where the latter integral vanishes exactly. For this reason In the next section we will show derivation of the corresponding system of equations, which holds, when the eigenfunctions are used.

2.1. On the spectral decomposition

In this section we compare the results of the previous section with the well-established ones which can be proved if the spectral decomposition, i.e., that made with eigenfunctions of the Stokes operator $\{\mathcal{W}_{k}\}$ is used, instead that of a generic POD basis. We recall that, by classical results about compact operators in Hilbert spaces there exists a sequence of smooth functions $\{\mathcal{W}_{k}\}$ (and their regularity is depending on the smoothness of the bounded domain $\Omega$) and an increasing sequence of positive numbers $\{\lambda_{k}\}$ such that

[TABLE]

Since each one the functions $\mathcal{W}_{k}$ solves the following Stokes system $A\mathcal{W}_{k}=\lambda_{k}\mathcal{W}_{k}$, then it follows by an integration by parts that

[TABLE]

hence also the $V$-orthogonality of the family $\{\mathcal{W}_{k}\}_{k\in\mathbbm{N}}$.

We consider now the usual decomposition by eigenfunctions associated with low and high frequencies

[TABLE]

where $\mathbf{P}_{m}$ is the projection over the subspace generated by $\{\mathcal{W}_{k}\}_{k=1,\dots,m}$. Our main result is based on a standard result about the projector $\mathbf{P}_{m}$, that can be found in [23, Appendix A.4, Theorem 4.11]:

Theorem 2.1**.**

The projector $\mathbf{P}_{m}$ can be defined as a continuous endomorphism over $V$, $H$ and $V^{\prime}$, and one has

[TABLE]

The result is mainly based on the regularity of solutions of elliptic equations, and thanks to this fact, it is possible to decompose the equations for the velocity, which yields,

[TABLE]

since $\mathbf{P}_{m}\Delta{\bf u}=\Delta\mathbf{P}_{m}\,{\bf u}=\Delta{\bf y}$ and also $\mathbf{Q}_{m}(\Delta{\bf u})=\Delta\mathbf{Q}_{m}\,{\bf u}=\Delta{\bf z}$.

Thus, we directly obtain the system

[TABLE]

which is more tractable than (2.2) from many points of view. In particular, the equations in (2.7), which can be rewritten also as

[TABLE]

do not contain the term $\mathcal{E}_{m}({\bf u})$.

3. Preliminaries on long-time averages

Since we consider long-time averages for the NSE, we must consider solutions which are global-in-time (defined for all positive times). Due to the well-known open problems related to the NSE, this forces us to restrict ourselves to Leray-Hopf weak solutions [6, 22]. By using a then natural setting we take the initial datum ${\bf u}(0)$ in $H$. The classical Leray-Hopf result of existence (but not uniqueness) of a global weak solution ${\bf u}$ to the NSE holds when ${\bf f}\in V^{\prime}$, and the velocity ${\bf u}$ satisfies

[TABLE]

Notice that consider in this paper the case where ${\bf f}$ is time-independent, for simplicity. However, the following results can be extended to the case where ${\bf f}={\bf f}(t)$ is time dependent, for ${\bf f}$ belonging to a suitable class (see [3]).

In order to properly set what we mean by long-time-averaging, let $\psi:\ \mathbbm{R}^{+}\times\Omega\to\mathbbm{R}^{N}$ be any tensor field related to a given turbulent flow ($N$ being its order). The time-average over a time interval $[0,t]$ is defined by

[TABLE]

According to the standard turbulence modeling process, we then apply the averaging operator $M_{t}$ to NSE (1.1) and also to (2.1), to study the limits when $t\to+\infty$. We recall that the long-times averages represent one of the few observable and computable quantities associated to a highly variable turbulent flow. We will adopt the following standard notation for the long-time average of any field $\psi$

[TABLE]

whenever the limit exists. (Without too much restrictions we can suppose that the limits we write do exist, at least after extracting sub-sequences leaving the mathematical difficulties, which can be treated with generalized Banach limits, for a more general and abstract framework). Within this theory we can decompose the velocity as follows

[TABLE]

where ${\bf u}^{\prime}$ represent the so called turbulent fluctuations. We recall that time-averaging has been introduced by O. Reynolds [29], at least for large values of $t$, and the ideas have been widely developed by L. Prandtl [24] in the case of turbulent channel flows. The same ideas have been also later considered in the case of fully developed homogeneous and isotropic turbulence, such as grid-generated turbulence. In this case the velocity field is postulated as oscillating around a mean smoother steady state, see for instance G.-K. Batchelor [1]. For further details on the role of time averaging in turbulence, after the work of Stokes and Reynolds, we can recall a few recent papers and books [2, 3, 5, 10, 16, 19, 21], where aspects of computation and modeling are studied.

We now observe that, by taking the time-averages of the NSE we have the following estimates, see [21, Prop. 2.1]

[TABLE]

where $C_{P}$ is the best constant in the Poincaré inequality

[TABLE]

The above inequalities show that both $\|{\bf u}(t)\|^{2}$ and $\frac{1}{t}\int_{0}^{t}\|\nabla{\bf u}(s)\|^{2}\,ds$ are uniformly bounded, hence we have the following result

Theorem 3.1** (cf. [3, 21]).**

Let ${\bf u}(0)\in H$, ${\bf f}\in V^{\prime}$, and let ${\bf u}$ be a global-in-time weak solution to the NSE (1.1). Then, there exist

1. (1)

a sequence $\{t_{n}\}_{n\in\mathbbm{N}}$ such that $\lim\limits_{n\to\infty}t_{n}=+\infty$; 2. (2)

a vector field $\overline{{\bf u}}\in V$; 3. (3)

a vector field ${\bf B}\in L^{3/2}(\Omega)^{3}$; 4. (4)

a second order tensor field $\boldsymbol{\sigma}^{(\textsc{r})}\in L^{3}(\Omega)^{9}$;

such that it holds:

1. i)

When $n\to\infty$,

[TABLE] 2. ii)

The Reynolds averaged equations:

[TABLE]

hold true in the weak sense; 3. iii)

The equality ${\bf F}={\bf B}-(\overline{{\bf u}}\cdot\nabla)\,\overline{{\bf u}}=\nabla\cdot\boldsymbol{\sigma}^{(\textsc{r})}$ is valid in $\mathcal{D}^{\prime}(\Omega)$; 4. iv)

The following energy balance (equality) holds true

[TABLE] 5. v)

The tensor $\boldsymbol{\sigma}^{(\textsc{r})}$ is dissipative in average or, more precisely, the following inequality

[TABLE]

holds true.

It is important to observe that the long-time limit is characterized by the solution of the system (3.3), which is similar to the Navier-Stokes equations, but which contains an extra term, coming from the effect of fluctuations, which has the mean effect of increasing the dissipation.

We observe that this is related to the long-time behavior of solutions close to statistical equilibrium. The study of the long-time behavior dates back to pioneering works of Foias and Prodi on deterministic statistical solutions, see for instance [11]. Their interest is mainly devoted to finding measure in the space of initial data to be connected with the long-time limits. Here, we follow a slightly different path, as in [3, 21], in order to characterize in a less technical way the long-time behavior, without resorting to any ergodic-type result and also with the perspective that long time averages are computable or at least can be approximated in a clear way.

4. Average transfer of energy at equilibrium

Our goal is now to characterize in some sense the energy transfer between the two functions ${\bf y}$, ${\bf z}$ of the expansion and to determine –if possible– the sign of $\mathrm{e}_{m}({\bf u})$, at least in an average sense.

The point concerning the exchange of energy between low and high modes is in the same spirit as the results recalled in Foias, Manley, Rosa, and Temam [11, Chap. 5] and follows from results obtained in a more heuristic way by Kolmogorov [17].

We first observe that the $L^{2}$-orthogonality of the POD decomposition implies that

[TABLE]

Hence, from the uniform $L^{2}$-bound on ${\bf u}$ it follows that both ${\bf y}$ and ${\bf z}$ are uniformly bounded in time. From this observation we can deduce the following result, reminding that $\mathrm{e}_{m}$ and and $\mathcal{E}_{m}$ are defined by equations (2.4) and (2.5), and $M_{t}$ is defined by equation (3.1).

Theorem 4.1**.**

There exist a sequence $\{t_{n}\}$ such that $t_{n}\to+\infty$ and a field $\overline{{\bf z}}\in H$ such that

[TABLE]

and

[TABLE]

Proof.

Let us observe first that by the energy inequality (3.2), we easily deduce that $(M_{t}({\bf z}))_{t>0}$ is bounded in $H$, hence the first assertion of the statement and (4.1). We next prove (4.2). To do so, we average with respect to time with the operator $M_{t}$ the balance equation (2.6) for $\mathrm{E}({\bf z})$, which yields

[TABLE]

By using the energy inequality (3.2) once again, we see that the first two terms vanish as $t\to+\infty$ and also that $M_{t}(\|\nabla{\bf z}\|^{2})$ is bounded. Therefore, (4.3) yields

[TABLE]

hence (4.2). We observe that in this case we do not have any direct estimation on the behavior of the $H^{1}$-norm. ∎

In the case we can assume that the limit exists, we also have the following result.

Corollary 4.1**.**

Let us assume the limit of $M_{t_{n}}(\mathrm{e}_{m}({\bf u}))$ for $n\to+\infty$ exists, and that

[TABLE]

Then, it follows

[TABLE]

This result can be interpreted as that, beyond the range of injection of energy, the average net transfer of energy occurs only into the small scales. This occurs if the term of interaction between gradients of large and small scales is negligible, in the limit of long times. This latter assumption is not proved rigorously, but we will see it is satisfied in the numerical tests, with a good degree of approximation (see Section 5). However, when one uses the eigen-vectors of the Stokes operator as POD basis, this is automatically satisfied since this basis is also orthogonal for the $H^{1}$-scalar product, so that in this case $\mathcal{E}_{m}({\bf u})=0$.

4.1. The spectral case

The results of the previous section can be made much more precise in the case of decomposition made by a spectral basis of eigenfunctions of the Stokes operator. We present the results, which are in some sense new and not fully completely included in [11], in the sense of time-averaging. This procedure is applied to ${\bf u}$, which is a weak solution of the Navier-Stokes equations, satisfying the uniform estimates (3.2). In this way, the orthogonality (in both $H$ and $V$) of the basis implies that

[TABLE]

The results in this case are more precise than those from Theorem 4.1, since we have at disposal a larger set of a priori estimates and also the set of equations (2.7) has a better structure than (2.2).

We now prove the following results in the case of a decomposition of the velocity into small and large frequencies. The first one aims at taking the time average and then let $t$ go to infinity in the equations (2.2) satisfied by ${\bf y}$ and ${\bf z}$. The second one aims at comparing the amount of turbulent dissipation of small and large frequencies with respect to the total work of the corresponding Reynolds stresses $\boldsymbol{\sigma}^{(\textsc{r})}_{\bf y}$ and $\boldsymbol{\sigma}^{(\textsc{r})}_{\bf z}$.

Theorem 4.2**.**

Let ${\bf u}(0)\in H$, ${\bf f}\in\mathbf{P}_{m}\,H$, and let ${\bf u}$ be a global-in-time weak solution to the NSE (1.1). Then, there exist

1. (1)

a sequence $\{t_{n}\}_{n\in\mathbbm{N}}$ such that $\lim\limits_{n\to\infty}t_{n}=+\infty$; 2. (2)

vector fields $\overline{{\bf y}},\overline{{\bf z}}\in V$; 3. (3)

vector fields ${\bf B}_{1},{\bf B}_{2}\in V^{\prime}$;

such that it holds:

1. i)

When $n\to\infty$,

[TABLE] 2. ii)

The Reynolds averaged equations:

[TABLE]

and

[TABLE]

holds true in $V^{\prime}$.

Arguing as in [3, 21], using the relations $(\overline{{\bf z}}\cdot\nabla)\,\overline{{\bf z}}=\nabla\cdot(\overline{{\bf z}}\otimes\overline{{\bf z}})$ and $(\overline{{\bf y}}\cdot\nabla)\,\overline{{\bf y}}=\nabla\cdot(\overline{{\bf y}}\otimes\overline{{\bf y}})$, we get the existence of “small frequencies” and “large frequencies” Reynolds stresses $\boldsymbol{\sigma}^{(\textsc{r})}_{\bf y}$ and $\boldsymbol{\sigma}^{(\textsc{r})}_{\bf z}$ in $V^{\prime}$, such that

[TABLE]

or equivalently, if we write the Reynolds decomposition as

[TABLE]

then

[TABLE]

where the bar operator denotes the limit of the $M_{t_{n}}$’s in $V^{\prime}$ as $n\to\infty$ (eventually after having extracted another sub-sequence).

According to the budget (3.4), we aim to compare the turbulent dissipation of small and large scales, denoted as $\varepsilon^{\downarrow}$ and $\varepsilon^{\uparrow}$ respectively, to the total work of the Reynolds stresses, namely $(\nabla\cdot\boldsymbol{\sigma}^{(\textsc{r})}_{\bf y},\overline{{\bf y}})$ and $(\nabla\cdot\boldsymbol{\sigma}^{(\textsc{r})}_{\bf z},\overline{{\bf z}})$, where

[TABLE]

When we compare to (3.4), we observe that triad nonlinear effect between small and large frequencies will be felt, that means the nonlinear interactions due to the convection will be provided by the term

[TABLE]

Notice that due to the regularity of ${\bf y}$, it is easily checked that the following energy balance holds true (this property will be shown with more details in the proof of Theorem 4.2 below)

[TABLE]

Finally, we will use the following orthogonality relation (see e.g. [3, Lemma 4.4]), formally written as follows

[TABLE]

Theorem 4.3**.**

The families $(M_{t}(\Phi_{{\bf z}}({\bf y})))_{t>0}$ and $(M_{t}(\Phi_{{\bf y}}({\bf z})))_{t>0}$ converge (along certain subsequences) as $t\to\infty$. Let $\overline{\Phi_{{\bf z}}({\bf y})}$ and $\overline{\Phi_{{\bf y}}({\bf z})}$ denote the corresponding limits. One has

[TABLE]

and the following dissipation balances hold true

[TABLE]

Remark 4.1**.**

Notice that by equations (4.8) and (4.9) we see that $\boldsymbol{\sigma}^{(\textsc{r})}_{\bf y}$ is dissipative in mean, and follows the same law (3.4) as the complete Reynolds stress, namely

[TABLE]

However, nothing similar can be concluded from (4.9) about $\boldsymbol{\sigma}^{(\textsc{r})}_{\bf z}$, that might be at this stage non dissipative in mean, which permits an inverse energy cascade to occur.

The results of Theorems 4.2 and 4.3 are original, even if similar results have been already obtained in [10] and reported also in [11]. In that case, the results are based on the notion of deterministic statistical solutions and on a sort of ergodic hypothesis. Even if statements could look very similar to ours, the main difference is that we do not average over the set $H$ of initial data, and we do not introduce probability measures on $H$, as suggested by the work by Prodi [25, 26]. Our approach is based on a more elementary functional setting and also amenable to include treatment of sets of external forces, as those in several numerical or practical experiments. The main point is an extension of the results in [3].

Proof of Theorem 4.2.

We know, from the results in [3, 21] that $\mathbf{U}_{t}=M_{t}({\bf u})$ is such that

[TABLE]

hence, if we define $\mathbf{F}:=\mathbf{B}-(\overline{{\bf u}}\cdot\nabla)\,\overline{{\bf u}}$, we get

[TABLE]

and using $\overline{{\bf u}}\in V$ as test function we obtain

[TABLE]

We assume now that $\mathbf{P}_{m}{\bf f}={\bf f}$, and we consider the equations satisfied by ${\bf y}=\mathbf{P}_{m}{\bf u}$ and ${\bf z}=\mathbf{Q_{m}}{\bf u}$. In particular, the equation for ${\bf y}$ reads, as an abstract equation in $V_{m}=\mathbf{P}_{m}V$, as follows:

[TABLE]

The uniform estimates on ${\bf u}$ from Theorem 3.1 combined with Theorem 2.1, about the properties of the projection operator $\mathbf{P}_{m}$ as a continuous operator over $V^{\prime}$, yield

[TABLE]

in such a way that $\overline{{\bf y}}$ satisfies

[TABLE]

where $\mathbf{F}_{{\bf y}}:=\mathbf{B}_{1}-(\overline{{\bf y}}\cdot\nabla)\,\overline{{\bf y}}$, which leads to (4.4) by De Rham Theorem. Arguing as in [3] (which was already mentioned above), it is easily checked that there exists $\boldsymbol{\sigma}^{(\textsc{r})}_{\bf y}$ such that $\mathbf{F}_{{\bf y}}=\nabla\cdot\boldsymbol{\sigma}^{(\textsc{r})}_{\bf y}$. Hence, being $\overline{{\bf y}}\in V$ a legitimate test function, we get

[TABLE]

The other term ${\bf z}$ of the decomposition satisfies

[TABLE]

The uniform estimates on ${\bf u}$ and the boundedness of $\mathbf{P}_{m}$ imply the following convergence (up to a sub-sequence), as already shown in Theorem 4.1

[TABLE]

By using that $\mathbf{B}=\mathbf{P}_{m}\mathbf{B}+\mathbf{Q}_{m}\,\mathbf{B}$, we get

[TABLE]

for $\mathbf{F}_{{\bf z}}=\mathbf{B}_{2}-(\overline{{\bf z}}\cdot\nabla)\,\overline{{\bf z}}=\nabla\cdot\boldsymbol{\sigma}^{(\textsc{r})}_{\bf z}$. Hence, (4.5) follows again by De Rham Theorem. Notice that, by using as test function $\overline{{\bf z}}\in V$ we have the following energy equality:

[TABLE]

which concludes the proof. ∎

Proof of Theorem 4.3.

We now write the energy inequality for ${\bf z}$, obtaining

[TABLE]

and hence, by using the orthogonality of the basis, we have that $\mathbf{Q}_{m}{\bf z}={\bf z}$ and

[TABLE]

recalling the definition of $\Phi_{{\bf z}}({\bf y})$ in (4.6).

Averaging over a fixed time interval $(0,t)$, we get

[TABLE]

The $L^{2}$-uniform bounds on ${\bf z}$ imply that $\frac{1}{2t}\|{\bf z}(t)\|^{2}\to 0$, hence, possibly after having extracted another sub-sequence to ensure the convergence of the term $M_{t}(\|\nabla{\bf z}\|^{2})$ (that is known to be bounded by the energy inequality (3.2)) we get

[TABLE]

We now combine the orthogonality relation (4.7) with the energy balance (4.12), so that (4.13) yields

[TABLE]

which is almost inequality (4.10), up to the convergence of $(M_{t}(\Phi_{{\bf z}}({\bf y})))_{t>0}$ that remains to be proved. To prove this, we deal with the budget for ${\bf y}$, reminding that

[TABLE]

Then, averaging the energy equality (in this case we have equality since ${\bf y}$ solves a finite dimensional system of ordinary differential equations) which is satisfied for ${\bf y}$, we get

[TABLE]

By the same argument, eventually after having extracted a further sub-sequence, $(M_{t}(\|\nabla{\bf y}\|^{2}))_{t>0}$ is convergent as $t_{n}\to\infty$, as well as $(<{\bf f},M_{t}({\bf y})>)_{t>0}$. Therefore, $\{M_{t_{n}}(\Phi_{{\bf y}}({\bf z}))\}$ is also convergent by (4.15). Let $\overline{\Phi_{{\bf y}}({\bf z})}$ denotes its limit. In particular, by (4.14), $\{M_{t_{n}}(\Phi_{{\bf z}}({\bf y})))\}$ is also convergent, with limit $\overline{\Phi_{{\bf z}}({\bf y})}=-\overline{\Phi_{{\bf y}}({\bf z})}$. We are done with (4.10).

It remains to check (4.9). Taking the limit as $n\to\infty$ in (4.15) gives the equality

[TABLE]

which, combined with the energy balance (4.11) and the orthogonality relation (4.7), yields (4.9), ending the proof. ∎

5. Numerical results

In Theorem 4.1, we showed that

[TABLE]

In this section, we investigate numerically whether the inequality (5.1) holds. To this end, we consider the one-dimensional Burgers equation with homogeneous boundary conditions as a simplified, yet relevant test case:

[TABLE]

To calculate the long-time average of $\mathrm{e}_{m}(u)$ in (5.1), we use the composite trapezoidal rule:

[TABLE]

where $t_{i}=(i-1)*\frac{T}{n},~{}i=1,...,n+1.$ We also use the composite trapezoidal rule to calculate the long-time average of $\mathcal{E}_{m}(u)$.

5.1. Numerical Results with step function initial condition

Our numerical results are obtained by using the one-dimensional Burgers equation (5.2) with a step function initial condition [15, 34]:

[TABLE]

We use the following parameters in the finite element discretization of the Burgers equation (5.2): $\Omega=[0,1]$, $\nu=10^{-2}$, $f=0$, mesh size $h=1/128$, piecewise linear finite element spatial discretization, and backward Euler time discretization.

5.1.1. Case 1:

For this test case, we consider the time interval $[0,T]=[0,1]$ and the time step $\Delta t=10^{-2}$. We utilize all the snapshots to build the POD basis, whose dimension is $d=37$. In the composite trapezoidal rule, we use $n=100$. In Figure 5.1.1, we plot the DNS results (which are used to generate the snapshots). In Table 1, we list the time-averages of $\mathrm{e}_{m}(u)$ and $\mathcal{E}_{m}(u)$ for different $m$ values. We note that the time-average of $\mathrm{e}_{m}(u)$ is positive for all $m$ values. The time-average of $\mathcal{E}_{m}(u)$ is positive for the low $m$ values and negative for the largest $m$ values. Furthermore, the magnitude of the time-average of $\mathcal{E}_{m}(u)$ is lower than the magnitude of the time-average of $\mathrm{e}_{m}(u)$. Thus, we conclude that the time average $\frac{1}{T}\int_{0}^{T}\mathrm{e}_{m}(u(s))+\mathcal{E}_{m}(u(s))\,ds$ in (5.1) is positive for all $m$ values.

In Case 1, we showed that the time average $\frac{1}{T}\int_{0}^{T}\mathrm{e}_{m}(u(s))+\mathcal{E}_{m}(u(s))\,ds$ in (5.1) is positive. In the remainder of this section, we investigate whether this time average remains positive if we make the following changes in our computational setting: (i) we increase/decrease the time-interval; and (ii) we use more quadrature points (i.e., subintervals) in the composite trapezoidal rule (5.3).

5.1.2. Case 2:

In this case, we use a longer time interval, i.e., $[0,T]=[0,10]$ (instead of $[0,T]=[0,1]$, as we used in Case 1). We also use different time step ($\Delta t$) values to generate the snapshots and a different number of quadrature points to evaluate the time average $\frac{1}{T}\int_{0}^{T}\mathrm{e}_{m}(u(s))$.

In Tables 3–5, we list the time-averages of $\mathrm{e}_{m}(u)$ and $\mathcal{E}_{m}(u)$ for different time steps ($\Delta t$) values, different number of equally spaced quadrature points ($n$), and different $m$ values. We note that the time-averages of $\mathrm{e}_{m}(u)$ and $\mathcal{E}_{m}(u)$ are positive for all $\Delta t$, $n$, and $m$ values. Furthermore, the magnitude of the time-average of $\mathcal{E}_{m}(u)$ is lower than the magnitude of the time-average of $\mathrm{e}_{m}(u)$. Thus, we conclude that the time average $\frac{1}{T}\int_{0}^{T}\mathrm{e}_{m}(u(s))+\mathcal{E}_{m}(u(s))\,ds$ in (5.1) is positive for all $\Delta t$, $n$, and $m$ values. Furthermore, we note that decreasing the time step while keeping the same number of snapshots (i.e., $n=10000$) does not change the time average $\frac{1}{T}\int_{0}^{T}\mathrm{e}_{m}(u(s))+\mathcal{E}_{m}(u(s))\,ds$ significantly (see Tables 3–5).

5.1.3. Case 3:

In this case, we use an even longer time interval, i.e., $[0,T]=[0,100]$, and compare the time-averages for this time interval to those for the time intervals $[0,T]=[0,1]$ (Case 1) and $[0,T]=[0,10]$ (Case 2). For each time interval, we use the same time step values ($\Delta t=10^{-2}$) to generate the snapshots and all the subintervals in the composite trapezoidal rule utilized in the evaluation of the time average $\frac{1}{T}\int_{0}^{T}\mathrm{e}_{m}(u(s))$. In Table 6, we list the time-averages of $\mathrm{e}_{m}(u)$ and $\mathcal{E}_{m}(u)$ for all three time intervals and different $m$ values. We note that the time-averages of $\mathrm{e}_{m}(u)$ and $\mathcal{E}_{m}(u)$ are positive for all time intervals and $m$ values. Furthermore, the magnitude of the time-average of $\mathcal{E}_{m}(u)$ is generally lower than the magnitude of the time-average of $\mathrm{e}_{m}(u)$. Thus, we conclude that the time-average $\frac{1}{T}\int_{0}^{T}\mathrm{e}_{m}(u(s))+\mathcal{E}_{m}(u(s))\,ds$ in (5.1) is positive for all time intervals and $m$ values. Furthermore, we note that the time-averages of $\mathrm{e}_{m}(u)$ and $\mathcal{E}_{m}(u)$ for the time intervals $[0,T]=[0,100]$ and $[0,T]=[0,10]$ are close, whereas those for the time interval $[0,T]=[0,1]$ are slightly different. Thus, we conclude that the time interval $[0,T]=[0,10]$ is adequate for the approximation of the long time-average $\frac{1}{T}\int_{0}^{T}\mathrm{e}_{m}(u(s))+\mathcal{E}_{m}(u(s))\,ds$.

5.1.4. Case 4:

In this case, we use a much shorter time interval, i.e., $[0,T]=[0,0.1]$, and compare the time-averages for this time interval to those for the time intervals $[0,T]=[0,1],[0,T]=[0,10],$ and $[0,T]=[0,100]$ (Case 3). We use two different time step values to generate the snapshots, but the same (i.e., $n=5000$) equally spaced subintervals in the composite trapezoidal rule utilized in the evaluation of the time average $\frac{1}{T}\int_{0}^{T}\mathrm{e}_{m}(u(s))$. In Tables 9–9, we list the time-averages of $\mathrm{e}_{m}(u)$ and $\mathcal{E}_{m}(u)$ for two different time step values and different $m$ values. We emphasize that, this time, the time-average of $\mathrm{e}_{m}(u)$ is negative for some $m$ values. Furthermore, the magnitude of the time-average of $\mathcal{E}_{m}(u)$ is this time larger than the magnitude of the time-average of $\mathrm{e}_{m}(u)$. This is in stark contrast with the previous cases.

6. Conclusions

In this preliminary study, we investigated theoretically and numerically the time-average of the exchange of energy among modes of reduced order models (ROMs) of fluid flows. In particular, we were interested in the statistical equilibrium problem, and especially in the long-time averaging of the ROM solutions. The main goal of the paper was to deduce the possible forward and backward average transfer of the energy among ROM basis functions (modes). We considered two types of ROM modes: eigenfunctions of the Stokes operator and proper orthogonal decomposition (POD) modes. In Theorem 4.1 and Theorem 4.2, we proved analytical results for both types of ROM modes and we highlighted the differences between them, especially those stemming from the lack of orthogonality of the gradients of the POD basis functions.

In Section 5, we investigated numerically whether the time-average energy exchange between POD modes (i.e., $\frac{1}{T}\int_{0}^{T}\mathrm{e}_{m}(u(s))+\mathcal{E}_{m}(u(s))\,ds$) in Theorem 4.1 is positive. To this end, we used the one-dimensional Burgers equation as a mathematical model. We utilized a piecewise linear FE spatial discretization and a backward Euler temporal discretization. To compute the time-averages, we used the composite trapezoidal rule. We tested different time steps, different number of subintervals in the composite trapezoidal rule, and, most importantly, different time intervals, to ensure that the computed quantities are indeed approximations of the time-averages and not numerical artifacts. The main conclusion of our numerical study is that, for long enough time intervals (i.e., time intervals longer than $[0,T]=[0,10]$), the time-average $\frac{1}{T}\int_{0}^{T}\mathrm{e}_{m}(u(s))+\mathcal{E}_{m}(u(s))\,ds$ in (5.1) is positive. Furthermore, the magnitude of the time-average of $\mathcal{E}_{m}(u)$ is much lower than the magnitude of the time-average of $\mathrm{e}_{m}(u)$.

There are several research directions that we plan to pursue. Probably the most important one is the numerical investigation of the theoretical results in three-dimensional, high Reynolds number flows, which could shed new light on the energy transfer among ROM modes. A related, but different numerical investigation was performed in [7].

Acknowledgments

Luigi C. Berselli and Roger Lewandowski for the research that led to the present paper were partially supported by a grant of the group GNAMPA of INdAM and by the project the University of Pisa within the grant PRA$\_{}2018\_{}52$ UNIPI Energy and regularity: New techniques for classical PDE problems. Traian Iliescu and Birgul Koc gratefully acknowledge the support provided by the NSF DMS-1821145 grant.