A mathematical model of stem cell regeneration with epigenetic state transitions

1.   Introduction

In this note we combine some results on the long-time averaging of fluid equations with the recent 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 and 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 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 averaging and link it 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)

with Dirichlet conditions when the motion takes place in a smooth and bounded domain $\Omega\subset \mathbb{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

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 the viscosity as a control parameter and 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, for limited computational and experimental resources, are the only ones computable and observable. 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 ^[2,5,27,39].

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 deep 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, to the computational practitioners. For these 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 much smaller approximate 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 ^{[15,25,35,36,38]}.

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

where $\{{\bf w}_k\}_{k\in \mathbb{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

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 generally used to determine the approximate velocity and also the energy content in the basis construction. To determine the number $r\in \mathbb{N}$ such that the solution is projected over the space 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 ^[41,42,44]. The present paper combines mathematical results on the long-time behavior of fluid flows (especially in the case of statistical equilibrium) with reduced order modeling. The main goal of this approach is that it provides a mathematical description of both the long-time averages of ROMs and the energy exchange between ROM modes. Furthermore, preliminary numerical results that support the theoretical developments are presented. Specifically, we are extending to the POD setting the results for statistical solutions by Foias et al. ^[11,12] and those more recently obtained for time-averages in ^[3,28]. 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,28]. These results show that low frequency modes yield a Reynolds stress that is dissipative in the mean, its total spatial mean work being larger than the long time average of the dissipation of the fluctuations, which is consistent with observations and results in ^[3,11]. 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.

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 perform a preliminary numerical study in which 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 a meaningful physical quantity and the analysis of its budget is at the basis of abstract existence results for weak solutions (cf. Constantin and Foias ^[6]) and also of the conventional turbulence theories of Kolmogorov ^[19]. 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

At present, it is possible to prove that the above balance holds true with the sign of "less or equal" for the class of weak (or turbulent) solutions for which there are results of existence 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

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

holds (at least with the inequality sign) for all smooth scalar functions $\phi\in C^\infty_0((0, T)\times\Omega)$ such that $\phi\geq0$. The validity of such an inequality is one of the requirements to prove partial regularity results, but it is also one of the conditions to be satisfied by weak solutions constructed by numerical or LES methods. In this respect, see Guermond, Oden and Prudhomme ^[14] and ^[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

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 ^[12], if we define

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

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:

(ⅰ) choose modes $\{ {\bf w}_1, \ldots, {\bf w}_d\}$, which represent the recurrent spatial structures of the given flow;

(ⅱ) choose the dominant modes $\{{\bf w}_1, \ldots, {\bf w}_m\}$, with $m \leq d$, as basis functions for the ROM;

(ⅲ) use a Galerkin truncation ${\bf u}_m = \sum_{j = 1}^{m} a_j \, {\bf w}_j$;

(ⅳ) replace ${\bf u}$ with ${\bf u}_m$ in the NSE;

(ⅴ) 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:

where $a$ is the vector of unknown ROM coefficients and $A, B$ are ROM operators;

(ⅵ) in an offline stage, compute the ROM operators;

(ⅶ) in an online 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 more energetic modes (${\bf y}$), and the remainder (${\bf z}$). This can be formalized as follows:

where

In (2.1), $m\in \mathbb{N}$ 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 in ${\bf y}$. Furthermore, ${P}_m$ is the projection operator over the subspace $V_m$ spanned by the first $m$-functions $\{{\bf w}_k\}_{k = 1, \dots, m}$.

We observe that we are considering the functions $\{{\bf w}_k\}_{k}$ as divergence-free, at least in a weak and/or approximate sense. Even if generally they are not "exactly divergence-free", numerically we can consider that they have vanishing divergence, by assuming that the solution of the problems used to construct the basis is accurate enough to have negligible divergence. Hence, in the computations that will follow, we will drop the pressure terms by a standard Leray projection. It will be nevertheless interesting to extend our study to bases that are not divergence-free, e.g., those derived by the combination of ROM with artificial compressibility methods ^[8,9,13].

In addition, we consider the external force as stationary, that is ${\bf f} = {\bf f}(\mathbf{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

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

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$

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 ^[12]), 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, in particular by Theorem 2.1), hence

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

where we used that both ${P}_m$ and ${Q}_m$ commute with the time derivative.

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

In this way we can obtain the following equality and inequality

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.3)

This is a formal setting, which is obviously true for strong solutions of the NSE, where the inequality in (2.4) 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.4) 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

We also use the following notation:

Hence, we can rewrite (2.4) as follows

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.2) and the long-time integration of its solution pose hard numerical problems.

We point out for the reader that 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 the 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 that can be proved if the spectral decomposition (i.e., that made with eigenfunctions of the Stokes operator $\{\mathcal{W}_k\}$) is used instead of utilizing 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

Since each function $\mathcal{W}_k$ solves the following Stokes system $A\mathcal{W}_k = \lambda_k \mathcal{W}_k$, it follows by an integration by parts that

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

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

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 [30,App. A.4,Thm. 4.11]:

Theorem 2.1. The projector $\mathbf{P}_m$ can be defined as a continuous endomorphism over $V$, $H$ and $V'$, and one has

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,

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

which can be rewritten also as

We note that the equations in (2.9) do not contain the term $\mathcal{E}_{m}({\bf u})$, whereas the equations in (2.7) contain $\mathcal{E}_{m}({\bf u})$. This will have important consequences in the analysis in Section 4.

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,29]. By using a 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'$, and the velocity ${\bf u}$ satisfies

Notice that we 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:\ \mathbb{R}^+\times\Omega\to \mathbb{R}^N$ be any tensor field related to a given turbulent flow ($N$ being its order). The time-average of $\psi$ over a time interval $[0, t]$ is defined by

According to the standard turbulence modeling process, we then apply the averaging operator $M_t$ to NSE (1.1) and also to (2.2), to study the limits when $t\to+\infty$. We recall that the long-time 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$

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

where ${\bf u}': = {\bf u}-\overline{{{\bf u}}}$ represents the so-called turbulent fluctuations. We recall that time-averaging has been introduced by O. Reynolds ^[37], at least for large values of $t$, and the ideas have been widely developed by L. Prandtl ^[32] 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,11,18,26,28]}, 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 [28,Prop. 2.1]

where $C_P$ is the best constant in the Poincaré inequality

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 for all $t\geq 1$ (any other positive time will be enough), hence we have the following result:

Theorem 3.1 (cf. ^[3,28]). Let ${\bf u}(0) \in H$, ${\bf f}\in V'$, and let ${\bf u}$ be a global-in-time weak solution to the NSE (1.1). Then, there exist

1. a sequence $\{t_n\}_{n\in \mathbb{N}}$ such that $\lim\limits_{n \to \infty }t_n = +\infty$;

2. a vector field $\overline{{{\bf u}}}\in V$;

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

4. a second order tensor field $\boldsymbol {\sigma}^{(\text{R})} \in L^3(\Omega)^9$;

such that it holds:

i) When $n \to \infty$,

ii) The Reynolds averaged equations:

hold true in the weak sense;

iii) The equality ${\bf B} - (\overline{{{\bf u}}} \cdot \nabla)\, \overline{{{\bf u}}} = \nabla \cdot \boldsymbol {\sigma}^{(\text{R})}$ is valid in $\mathcal{D}'(\Omega)$;

iv) The following energy balance (equality) holds true

v) The tensor $\boldsymbol {\sigma}^{(\text{R})}$ is dissipative on the average or, more precisely, the following inequality

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 ^[12]. 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,28], 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 in this section is to characterize in some sense the energy transfer between the two functions ${\bf y}$ and ${\bf z}$ of the expansion and to determine –if possible– the sign of $\mathrm{e}_m({\bf u})$, at least in an average sense. We will consider both the POD case and the spectral one.

4.1.   The POD case

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 [12,Chap. 5] and follows from results obtained in a more heuristic way by Kolmogorov ^[19].

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

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 $\mathcal{E}_{m}$ are defined by equations (2.5) and (2.5), and $M_t$ is defined by equation (3.1).

Theorem 4.1. Let ${\bf z}$ be the projection onto the less energetic POD modes of a weak solution ${\bf u}$ to the Navier-Stokes equations. Then, there exists a sequence $\{t_n\}$ such that $t_n\to+\infty$ and a field $\overline{{{\bf z}}}\in H$ such that

and

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.7) for $\mathrm{E}({\bf z})$, which yields

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

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

Then, it follows

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 time. 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 eigenfunctions 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.2.   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 ^[12], 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

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.8) has a better structure than (2.3).

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.3) 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}^{(\text{R})}_{\bf y}$ and $\boldsymbol {\sigma}^{(\text{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. a sequence $\{t_n\}_{n\in \mathbb{N}}$ such that $\lim\limits_{n \to \infty }t_n = +\infty$;

2. vector fields $\overline{{{\bf y}}}, \overline{{{\bf z}}}\in V$;

3. vector fields ${\bf B}_1, {\bf B}_2\in V'$;

such that it holds:

i) When $n \to \infty$,

ii) The Reynolds averaged equations:

and

holds true in $V'$.

Arguing as in ^[3,28], 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}^{(\text{R})}_{\bf y}$ and $\boldsymbol {\sigma}^{(\text{R})}_{\bf z}$ in $V'$, such that

or equivalently, if we write the Reynolds decomposition as

then

where the bar operator denotes the limit of the $M_{t_n}$'s in $V'$ 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}^{(\text{R})}_{\bf y}, \overline{{{\bf y}}})$ and $(\nabla \cdot\boldsymbol {\sigma}^{(\text{R})}_{\bf z}, \overline{{{\bf z}}})$, where

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

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)

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

which is valid for any field $\psi:\, \mathbb{R}^+\to V$, such that $\overline{{\psi}}$ is well-defined and the fluctuations are defined as $\psi' = \psi-\overline{{\psi}}$.

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

and the following dissipation balances hold true

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

However, nothing similar can be concluded from (4.9) about $\boldsymbol {\sigma}^{(\text{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 ^[11] and reported also in ^[12]. 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 ^[33,34]. 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,28] that $\mathbf{U}_t = M_t ({\bf u})$ is such that

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

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

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:

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'$, yield

in such a way that $\overline{{{\bf y}}}$ satisfies, for the spectral basis $\mathcal{W}_k$ introduced in Section 2.1,

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, if written in a strong formulation. Arguing as in ^[3] (which was already mentioned above), it is easily checked that there exists $\boldsymbol {\sigma}^{(\text{R})}_{\bf y}$ such that $\mathbf{F}_{{\bf y}} = \nabla \cdot \boldsymbol {\sigma}^{(\text{R})}_{\bf y}$. Hence, being $\overline{{\bf y}}\in V_m\subset V$ a legitimate test function, we get

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

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

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

for $\mathbf{F}_{{\bf z}} = \mathbf{B}_2-(\overline{{\bf z}}\cdot\nabla)\, \overline{{\bf z}} = \nabla \cdot \boldsymbol {\sigma}^{(\text{R})}_{\bf z}$. Hence, (4.5) follows again by De Rham Theorem. Notice that, $\overline{{{\bf z}}}\in V_m^\perp\subset V$ and in particular $\overline{{{\bf z}}}\in \text{span}\{\mathcal{W}_j\}_{j\geq m+1}$, so by linearity it can be used as test function in (4.12). Next, since $\nabla\cdot\overline{{{\bf z}}} = 0$, it follows that $< (\overline{{\bf z}}\cdot\nabla)\, \overline{{\bf z}}, \overline{{{\bf z}}} > = 0$ and we have the following energy equality:

which concludes the proof.

Proof of Theorem 4.3. We now write the energy inequality for ${\bf z}$, obtaining

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

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

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

The $L^2$-uniform bounds on ${\bf z}$ imply that $\frac{1}{2t}\|{\bf z}(t)\|^2\to0$, 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

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

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}$, recall the definition of $\Phi_{{\bf z}}({\bf y})$ and $\Phi_{{\bf y}}({\bf z})$ from (4.6).

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

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})}}$ denote its limit. In particular, by (4.6), $\{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

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

In this section, we investigate numerically whether the inequality (5.1) holds. To this end, we consider the one-dimensional Burgers equation with homogeneous Dirichlet boundary conditions as a simplified, yet relevant test case ^{[16,17,20,21,22,23,24,31,40,43]}:

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

where $t_i = (i-1)*\frac{T}{n}, \; i = 1, \dots, 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 ^[16,43]:

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 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 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 $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: (ⅰ) we increase/decrease the time-interval; and (ⅱ) 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)) + \mathcal{E}_{m}(u(s)) \, ds$.

In Tables 2–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. Moreover, 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 $\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)) + \mathcal{E}_{m}(u(s)) \, ds$. 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)) + \mathcal{E}_{m}(u(s)) \, ds$. In Tables 8–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 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 performed a preliminary numerical investigation aiming at determining 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. We thank the two reviewers for the comments and suggestions that improved the paper.

Conflict of interest

The authors declare no conflict of interest.