On the Reynolds time-averaged equations and the long-time behavior of Leray-Hopf weak solutions, with applications to ensemble averages
Luigi C. Berselli Roger Lewandowski
math.AP







On the Reynolds time-averaged equations and the long-time behavior of
  Leray-Hopf weak solutions, with applications to ensemble averages
    Luigi C. Berselli
  

  Università di Pisa
  

  Dipartimento di Matematica
  

  Via Buonarroti 1/c, I-56127 Pisa, ITALY.
  

  email: luigi.carlo.berselli@unipi.it
  

  

  Roger Lewandowski
  

University of Rennes & INRIA
  

CNRS, IRMAR, UMR 6625 & Fluminance team, 
  

  Campus Beaulieu, F-35042 Rennes
  

  email: Roger.Lewandowski@univ-rennes1.fr 
    
==========================================================================================================================================================================================================================================================================================================================================================================


  We consider the three dimensional incompressible Navier-Stokes equations with non
  stationary source terms  chosen in a suitable space. We prove the existence of
  Leray-Hopf weak solutions and that it is possible to characterize (up to sub-sequences)
  their long-time averages, which satisfy the Reynolds averaged equations, involving a
  Reynolds stress. Moreover, we show that the turbulent dissipation is bounded by the sum
  of the Reynolds stress work and of the external turbulent fluxes, without any additional
  assumption, than that of dealing with Leray-Hopf weak solutions.

  Finally, in the last section we consider ensemble averages of solutions, associated to a
  set of different forces and we prove that the fluctuations continue to have a
  dissipative effect on the mean flow.



Keywords: Navier-Stokes equations, time-averaging, Reynolds equations,
Boussinesq hypothesis.



2010 MSC: Primary 35Q30. Secondary: 76D05, 76D06, 76F05.





§ INTRODUCTION


Let us consider the 3D homogeneous incompressible Navier-Stokes equations (NSE in the
sequel),

    {_t+ (·∇) -νΔ+∇
          p   =        in ]0,+∞[ ×Ω,
          
    ∇·   =0        in ]0,+∞[×Ω,
          
       =0        on ]0,+∞[×Γ,
          
    (0, )   =_0()         in Ω,
      .
 
where ⊂^3 is a bounded Lipschitz domain, = its boundary, =
(t, ) denotes the fluid velocity, p = p(t, ) the pressure, =  (t, )
the external source term, and (t, ) ∈ ]0,+∞[×Ω. The main aim of this
paper is to study the long-time averages of weak solutions to the NSE (<ref>),
namely 

    () := lim_t →∞ M_t (),   where   M_t() :=
      1/t∫_0^t  (s, )  ds,

when the source term  is time dependent, and to link it to the Reynolds averaged
equations (see (<ref>) below). We also will consider the problem of
ensemble averages, which is closely related to long-time averages.

Long-time average for “tumultuous” flows, today turbulent flows, seem to have been
considered first by G. Stokes <cit.> and then by O. Reynolds <cit.>. The

idea is that steady-state turbulent flows are oscillating around a stationary
flow, which can be expressed through long-time averages. L. Prandtl used them to introduce
the legendary “Prandtl mixing length” (see <cit.> and in <cit.>) to model a turbulent boundary layer over a plate.  However, although
long-time averaging plays a central role in turbulence modeling (since it is natural as
well as “ensemble averaging”) it is not clear whether the limit (<ref>) is
well defined, or not.

The mathematical problem of properly defining long-time averages and to investigate the
connection with the Reynolds equations was already studied before, but when the source
term  does not depend on time, namely =  ().

So far as we know, the first who considered this issue is
C. Foias <cit.>, when ∈ H, where


    H := {∈ L^2()^3  s.t.  · = 0, ·𝐧=0  on Γ}.
 
Foias analysis is based on the notion of “statistical solution” introduced
in <cit.> and he was able to prove that for any given Leray-Hopf
solution = (t,) to (<ref>) then:
 


- There exists ∈ H, such that, up to a sub-sequence, M_t () → in H as t →∞;

- There is a stationary statistical solution μ to the NSE, which is a probability
measure on H, such that it holds = ∫_H   dμ ();
 
- There exists a random variable on (H, μ) called ' (in the notation of
Foias called δ) such that

    '=∫_H' dμ=0    and    ∇'^2=∫_H∇'^2 dμ<∞,

such that (, ') is a solution to the Reynolds Equations given
by (<ref>) below. If μ is a statistical solution, ' is the random



variable expressed by its probability law, 

    Prob ( '∈ F) = μ (+ F),

for any Borel
set F∈ H.
 
 - According to our notations, the Reynolds stress  given by
 
    · = ∫_H ·[(-) ⊗ (-)] dμ (),

is dissipative on the mean flow, which is one of the main challenges of such
analysis, because of the Boussinesq assumption (see
in <cit.>). In fact this is much more better, since Foias
in <cit.> was able to prove that the turbulent
dissipation  is bounded by the work of the Reynolds stress on the mean flow,


    := ν ' ^2≤ (·, ),


where (·, ·) and . denote the standard L^2(Ω) scalar product
and norm, respectively. (Sometimes norm is normalized by |Ω|, the measure of the
domain, but this is inessential). This analysis is partially reported in <cit.>, where the limit in (<ref>) is replaced by the abstract
Banach limit, but without any link to the Reynolds equations and dissipation inequality
such as (<ref>). 

The original Foias analysis is very deep and essential to the field. However, it is worth
noting that in this approach:



i) The natural time filter used to determine the Reynolds stress (initially
suggested by Prandtl, when the limit makes sense) 


    =  () := lim_t →∞1/t∫_0^t [( - ) ⊗ ( - )](s, )  ds,


is replaced by an abstract probability measure that it is not possible to calculate in
practical simulations, although the ergodic assumption -which remains to be proved- would
mean that they coincide (see for instance Frisch <cit.>);



ii) The fluctuation given by (<ref>), when the force is
time independent, is a time independent random variable, which may be questionable from
the physical point of view. In fact, when  is given by (<ref>), one cannot
conclude that this ' yields the Reynolds decomposition

    (t, ) =  () + '(t, ),

in which the fluctuation is time dependent for a realistic non stationary flow.



iii) Concerning items i) and ii) above, we can suggest that probably there are
still work to be done in the interesting field of statistical solutions, if one wants to
use this mathematical tool to face the tough closure problem in the Reynolds equations,
and apply to realistic problem. Moreover, the case of a time dependent source term also
remains to be considered in connection with the scope of statistical solutions (see the
warning in <cit.>).
 


Still in the case of a stationary source term, the long-time averaged problem has been
also considered more recently in <cit.> and at the time the author was not aware
yet of the connection between time-averaging and statistical solutions in Foias
work. In <cit.>, he studied the equation satisfies by M_t() (see
equations (<ref>) below) for a Leray-Hopf weak solution of the NSE,
and he took the limit in this equation when the domain if of class C^9/4, 1 and
=()∈ L^5/4()^3 ∩ V', where


    V:= {∈ H^1_0()^3   s.t.  · = 0 },


and V' denotes the dual space of V, with duality pairing <. , .>.  The analysis is
based on the energy inequality, which yields a uniform estimate in time of the L^2 norm
of , and on a L^p-regularity result by Amrouche and Girault <cit.> about
the steady Stokes problems, valid in C^k, α domains. It is shown
in <cit.> that there exists ∈ L^5/3 ()^9 and ∈ W^1, 5/4
()/ such that, up to a sub-sequence, when t →∞, M_t () converges to a
field ∈ W^2, 5/4()^3, which satisfies in the sense of the distributions
the closed Reynolds equations:
 
    {
         (·∇)  -νΔ+∇p+∇·   =         in Ω,
          
    ∇·   =0        in Ω,
          
       =0        on Γ.
        .

Moreover, it is also shown that  is dissipative on the mean flow, namely
 
    0 ≤ (∇·, ),

which is weaker than (<ref>). 
 
The main part of the present study is in continuation of <cit.>, bringing
substantial improvements. The novelty is that we are considering time dependent source
terms = (t, ), which was never considered before for this problem, up to our
knowledge. Moreover, we do not need extra regularity assumption on the domain ⊂^3 and on . The first main result of this paper, Theorem <ref> below, is
close to that proved in <cit.>. Roughly speaking, we will show that M_t ()
given by (<ref>), converges to some  (up to sub-sequences) and that there
are  and  such that (<ref>) holds, at least in D'(), in which  is replaced by .
 
One key point is the determination of a suitable class for the source term. Throughout the
paper, we will take : _+→ V', made of function for which there is a
constant C>0 such that 

    ∀  t ∈_+    ∫_t^t+1 (s) ^2_V'  ds ≤ C.




In this respect we observe our results improve the previous ones also in terms of
regularity of the force, not only because we consider a time-dependent one. The main
building block of our work is the derivation of a uniform estimate of the L^2 norm of
, for  as above (see (<ref>) and its
Corollary (<ref>)). This allows to prove the existence of a weak
solution on [0,∞) to the NSE and to pass to the limit in the equation satisfied by
M_t(), when t →∞.

Moreover, we are able to generalise Foias result (<ref>) by proving
that the turbulent dissipation  is bounded by the sum of the work of the Reynolds
stress on the mean flow and of the external turbulent fluxes, namely


    = ν ' ≤ (·, ) + < ', ' > ,


which is one of the main features of our result.  Note that physically, it is expected
that (<ref>) becomes an equality for strong solutions.  Furthermore,
we also prove that when  is “attracted” (in some sense,
see (<ref>) below) by a stationary force = () as t →∞, then the turbulent fluxes < ', ' >
in (<ref>) vanish, so that (<ref>) is
restored. In particular, the question whether the stress tensor  is dissipative
remains an open problem for general unsteady , as those for which we still have global
existence of weak solutions.

In the second part of the paper, we will also consider ensemble averages, often used in
practical experiments. This consists of considering n∈ realizations of the flow,
{_1,…,_n} and to evaluate the arithmetic mean

    := 1/n∑_k=1^n _k.

Layton considers such ensemble averages in <cit.>, by introducing the
corresponding Reynolds stress, written as


    R(, ) =  ⊗ -⊗.


He shows that for a fixed stationary source term  and n strong solutions of the
NSE, then R(, ) is dissipative on the ensemble average, in time average. More
specifically it holds 


    lim inf_t →∞ M_t [ (· R(, ), ) ]  ≥ 0,


and to prove this inequality he first performs the ensemble average, then it takes
the long-time average.

In order to remove the additional assumption used in <cit.> of having strong
solutions, we will carry out an approach very close, but in a reversed order: We first take
the long-time average of the realizations for a sequence of time independent source terms
k. Then, we form the ensemble average

    𝐒^n:= 1/n∑_k=1^n _k,

where _k are weak solutions of the Reynolds equations.

Under suitable (but very light) regularity assumption about the _k's, we show the
convergence of S n to some <  > that satisfies the closed
Reynolds equations, and such that dissipativity still holds, that is


    0 ≤ ( · ,   ),
 
which holds for weak solutions (see the specific statement in
Theorem <ref>).



Plan of the paper.
The paper is organized as follows. We start by giving in Sec. <ref> the specific
technical statements of the results we prove. Then, in Sec. <ref> we give
some results about the functional spaces we are working with and we prove in detail the
main energy estimates (<ref>) and (<ref>) for source terms ∈ L^2_uloc(_+, V'), which yields an existence result of global weak solutions to
the NSE for such .  The Sec. <ref> is devoted to the Reynolds
problem and the develop additional properties of the time average operator M_t. Finally,
we give the proofs of the two main results on time averages in Sec. <ref>
and Sec. <ref>, respectively. 




§ MAIN RESULTS




 §.§ On the source term and an existence result
 

Since we aim to 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 enforces us to restrict to weak solutions. By using a most
natural setting, we take the initial datum _0 ∈ H, where H is defined
by (<ref>).



The classical Leray-Hopf results of existence (but without uniqueness) of a global weak
solution  to the NSE holds when ∈ L^2(_+;V'), and the velocity 
satisfies


    ∈ L^2(_+, V) ∩ L^∞(_+, H),


where V is defined by (<ref>) 
and V' denotes its topological dual. We will also denote by <  ,   > the duality
pairing[Generally speaking and when no risk of confusion occurs, we always denote
  by <  ,   > the duality pairing between any Banach space X and its dual X',
  without mentioning explicitly which spaces are involved.]  between V' and V.


Source terms ∈ L^2(_+;V') verify ∫_t^∞ f (s) _V'^2  ds
→ 0 when t →+∞. Therefore, a turbulent motion cannot be maintained for large
t, which is not relevant for our purpose. The choice adopted in the previous studies on
Reynolds equations was that of a constant force, and we also observe that many estimates
could have been easily extended to a uniformly bounded ∈ L^∞(_+;V'. On
the other hand, we consider a broader class for the source terms. According to the usual
folklore in mathematical analysis, we decided to consider the space L^2_uloc (_+;V')
made of all strongly measurable vector fields f: _+→ V' such that

    _L^2_uloc(_+;V')
       :=[sup_t≥0∫_t^t+1(s)_V'^2  ds ]^1/2<+∞.

 We will see in the following, that the above space, which strictly contains both
 L^2(_+;V') and L^∞(_+;V') is well suited for our framework. We will prove
 the following existence result, in order to make the paper self-contained.


   
  Let _0 ∈ H, and let f∈ L^2_uloc(_+;V'). Then, there exists a weak
  solution  to the NSE (<ref>) global-in-time, obtained by Galerkin
  approximations, such that
 
 
    ∈ L^2_loc (_+;V) ∩ L^∞(_+;H),


 and which satisfies for all t ≥ 0,
 
    (t) ^2 ≤_0 ^2 +  ( 3 + C_/ν )
        F^2/ν, 
       
       ν∫_0^t  (s) ^2 ds ≤_0 ^2 + ([t] + 1 ) 
         F^2/ν,

where ℱ:=_L^2_uloc(_+;V').


  The weak solution  shares most of the properties of the Leray-Hopf weak solutions,
  with estimates valid for all positive times.  Notice that we do not know whether or not
  this solution is unique. Anyway, it will not get regular as t →+∞, which
  is the feature of interest for our study.  As usual by regular we mean that it does not
  necessarily have the L^2-norm of the gradient (locally) bounded, hence that it is not
  a strong solution.





 §.§ Long-time averaging


This section is devoted to state the main results of the paper about long-time and
ensemble averages. 

  
  Let be given _0 ∈ H, ∈ L^2_uloc(_+;V'), and let  a
  global-in-time weak solution to the NSE (<ref>).  Then, there exist


  * a sequence {t_n}_n∈ such that lim_n →∞t_n=+∞;

  * a vector field ∈ V;

  * vector field ∈ V';

  * a vector field ∈ L^3/2(Ω)^3;

  * a second order tensor field ∈ L^3(Ω)^9;
 
 such that it holds:



  * when n →∞,
    
    M_t_n()⇀    in V,
            
        M_t_n () ⇀    in V',
            
        M_t_n ((·∇)  )⇀    in  L^3/2(Ω)^3,
            
       M_t_n('⊗' )⇀    in L^3(Ω)^9;
        M_t_n ( < ,  > ) →   <,  > +  < ', ' > ,

    where ' = -, and ' =  -;
  
  * the closed Reynolds equations (<ref>) holds true in the weak sense: 
  
    {
         (·∇)  -νΔ+∇p+∇·   =        in Ω,
          
    ∇·   =0        in Ω,
          
       =0        on Γ;
        .


  
  * the following equalities =  - (· )  = · are valid in 𝒟'(Ω);

  
  * the following energy balance holds true
  
    ν∇^2+∫_Ω· d=   <, >;



  * the turbulent dissipation  is bounded by the sum of the work of the Reynolds
  stress on the mean flow and the external turbulent fluxes,

    =  ν∇'^2≤∫_Ω(∇·)· d + < ', ' > ;



  *   if in addition the source term  verifies:


    ∃ ∈ V',   such that  lim_t→+∞∫_t^t+1(s)-_V'^2 ds=0,


then = and < ', ' >  = 0; in
particular, the Reynolds stress  is dissipative in average, that
is (<ref>) holds true.





Our second result has to be compared with results in Layton et al. <cit.>, where the long-time averages are taken for an ensemble of solutions.


  
  Let be given a sequence k⊂ L^q(Ω) converging weakly to some
   in L^q(Ω), with q>6/5 and let {^ k}_k ∈ be the associated long-time average of velocities, whose existence has been
  proved in Theorem <ref>. Then, the sequence of arithmetic averages of
  the long-time limits {^n}_n∈, defined as

    ^n :=  1/n∑_k=1^n^k,

converges weakly, as n→+∞, in V to some , which satisfies the
following system of Reynolds type

    {(·∇)  -νΔ+∇p+∇·   =        in Ω,
          
    ∇·   =0        in Ω,
          
       =0        on Γ,
      .

where  is dissipative in average, that is more precisely

    0 ≤1/||∫_Ω(∇·)· d.



In this case we do not have a sharp lower bound on the dissipation as in
Theorem <ref>, since here the averaging is completely different and the
fluctuations are not those emerging in long-time averaging. Nevertheless, the main
statement is in the same spirit of the first proved result.



§ NAVIER-STOKES EQUATIONS WITH UNIFORMLY-LOCAL SOURCE TERMS



This section is devoted to sketch a proof of Theorem <ref>.  Most
of the arguments are quite standard and we will give appropriate references at each step,
to focus on what seems (at least to us) non-standard when f∈
L^2_uloc(_+;V'); especially the proof of the uniform
L^2-estimate (<ref>), which is the building block for the results of the
present paper. Before doing this, we introduce the function spaces we will use, and
precisely define the notion of weak solutions we will deal with.



 §.§ Functional setting



Let Ω⊂^3 be a bounded open set with Lipschitz boundary ∂Ω. This is a sort of minimal assumption of regularity on the domain, in order to
have the usual properties for Sobolev spaces and to characterize in a proper way
divergence-free vector fields in the context of Sobolev spaces, see for instance
Constantin and Foias <cit.>, Galdi <cit.>, Girault and
Raviart <cit.>, Tartar <cit.>.

We use the customary Lebesgue spaces (L^p(Ω), . _p) and Sobolev spaces
(W^1,p(Ω), . _1,p).  For simplicity, we denote the L^2-norm simply by
. and we write H^1(Ω):=W^1,2(Ω). For a given sequence
{x_n}_n∈⊂ X, where (X, . _X) is Banach space, we denote by
x_n→ x the strong convergence, while by x_n⇀ x the weak one.

As usual in mathematical fluid dynamics, we use the following spaces


    V = {φ∈ D()^3,   ·φ = 0 }, 
      
       H =
        {∈ L^2(Ω)^3,  ∇·=0,·𝐧=0  on Γ},
        
       V = 
        {∈ H^1_0(Ω)^3,  ∇·=0},


and we recall that V is dense in H and V for their respective
topologies <cit.>.

Let (X, . _X) be a Banach space, we use the Bochner spaces L^p(I;X), for I=
[0,T] (for some  T>0) or I=_+ equipped with the norm


    u _L^p(I; X) := {   ( ∫_I  u(s) ^p_X ds  )^1/p     for 1≤ p<∞,
          
       ess sup_s∈ I u(s) _X     for p=+∞.
        .

 
The existence of weak solutions for the NSE (<ref>) is generally proved in the
literature when _0 ∈ H and the source term f∈ L^2(I; V'), or
alternatively when the source term is a given constant element of V'.  In order to study
the long-time behavior of weak solutions of the NSE (<ref>), we aim to enlarge the
class of function spaces allowed for the source term , to catch a more complex
behavior than that coming from constant external forces, as initially developed
in <cit.>. To do so, we deal with “uniformly-local” spaces, as defined below in
the most general setting.


  Let be given p∈[1,+∞[. We define L^p_uloc(_+;X) as the space of
  strongly measurable functions f:_+→ X such that


    f_L^p_uloc(X) :=[sup_t≥0∫_t^t+1f(s)_X^p ds ]^1/p<+∞.


It is easily checked that the spaces L^p_uloc(_+;X) are Banach spaces strictly containing
both the constant X-valued functions, and also L^p(_+; X), as illustrated
by the following elementary lemma.

  Let be given f∈ C(_+; X) converging to a limit ℓ∈ X, when
  t→+∞. Then, for any p∈[1,+∞[, we have that f ∈
  L^p_uloc(_+;X), and there exists T>0 such that



    f_L^p_uloc(X)≤ [sup_ t ∈ [0, T+1]f (t)_X^p + 2^p-1(1+ℓ_X)^p  ]^1/p .



 
  As ℓ = lim_t →+∞ f(t), there exists T >0 such that:
  ∀ t>T, f(t) -ℓ_X ≤ 1.  In particular, it holds


    ∫_t^t+1 f(s) ^p_X  ds ≤ 2^p-1(1+ ℓ_X)^p    for t>T,


while for all t ∈[0, T], 


    ∫_t^t+1 f(s) ^p_X ds  ≤sup_t ∈ [0, T+1]  f(t) _X^p,

hence the result.
 
However, it easy to find examples of discontinuous functions in L^p_uloc(_+;X)
which are not converging when t →+∞, and which are not belonging to
L^p(_+;X).



 §.§ Weak solutions


There are many ways of defining weak solutions to the NSE (see also
P.-L. Lions <cit.>). Since we are considering the incompressible case, the pressure
is treated as a Lagrange multiplier. Following the pioneering idea developed by
J. Leray <cit.>, the NSE are projected over spaces of divergence-free
functions. This is why when we talk about weak solutions the NSE, only the velocity 
is mentioned, not the pressure.

As in J.-L. Lions <cit.>, we give the following definition of weak solution, see
also Temam <cit.>.


  Given _0∈ H and f∈ L^2(I; V') we say that  is a weak solution
  over the interval I = [0,T] if the following items are fulfilled:
  
  
  * the vector field  has the following regularity properties

    
    ∈ L^2 (I; V) ∩ L^∞(I; H),


    and is weakly continuous from I to H, while
    lim_t→0^+(t)-_0_H=0;
  
  * for all φ∈ V,
  
    
    d/dt∫_ (t, ) ·φ()  d   - ∫_(t, ) ⊗(t, ) : φ() d
        + ν∫_ (t, ) : φ() d = 
            < f(t), φ > ,

    holds true in 𝒟'(I);
  
  * the energy inequality
    
    d/dt1/2(t)^2+ν∇(t)^2≤ <(t),(t) > ,

    holds in D'(I), where we write (t) instead of (t, ·) for
    simplicity.
  
  When ∈ L^2(0,T;V') and  is a weak solution in I=[0,T], and this holds
  true for all T>0, we speak of a global-in-time solution, or simply a global
    solution. In particular, ii) is satisfied in the sense of  𝒟'(0,+∞).


There are several ways to prove the existence of (at least) a weak solution to the NSE.
Among them, in what follows, we will use the Faedo-Galerkin method. Roughly speaking, let
{φ_n}_n∈ denote a Hilbert basis of V, and let, for
n∈, V_n := span{φ_1,⋯,φ_n}. By
assuming f∈ L^2(I; V'), it can be proved by the Cauchy-Lipschitz theorem
(see <cit.>) the existence of a unique _n ∈ C^1(I; V_n) such that for all
φ_k, with k= 1, …, n it holds


    d/dt∫__n (t, ) ·φ_k() d   -
        ∫__n(t, ) ⊗_n(t, ) : φ_k() d
        + ν∫__n (t, ) : φ_k() d =  <
          f(t), φ_k >,


 and which naturally satisfies the energy balance (equality)

 
    d/dt1/2_n(t)^2+ν∇_n(t)^2= <,_n >.


It can be also proved (always see again <cit.>) that from the sequence n one can extract a sub-sequence converging, in an appropriate sense, to a weak solution
to the NSE.  When I = _+ we get a global solution.

However, if assume f∈ L^2_uloc(_+;V'), the global result of existence does
not work so straightforward. Of course, for any given T>0, we have


    L^2_uloc(_+; V')_|[0,T]↪  L^2([0,T]; V'),


where L^2_uloc(_+; V')_|[0,T] denotes the restriction of a function in
L^2_uloc(_+; V') to [0,T]. Therefore, no doubt that the construction above holds
over any time-interval [0,T]. In such case letting T go to +∞ to get a global
solution (with some uniform control of the kinetic energy) is not obvious, and we do not
know any reference explicitly dealing with this issue, which deserves to be investigated
more carefully.  This is the aim of the next subsection, where we prove the most relevant
a-priori estimates.



 §.§ A priori estimates


Let be given f∈ L^2_uloc(_+; V'), and let _n = be the solution of
the Galerkin projection of the NSE over the finite dimensional space V_n. The function
 satisfies (<ref>) and (<ref>) (we do not write the subscript
n∈ for simplicity), is smooth, unique, and can be constructed by the Cauchy-Lipschitz
principle over any finite time interval [0,T]. Hence, we observe that by uniqueness it
can be extended to _+. We then denote F :=  f_ L^2_uloc (_+;
  V') and then after a delicate manipulation of the energy balance combined with the
Poincaré inequality, we get the following lemma.



For all t ≥ 0 we have

    (t)^2≤_0^2+(3+C_/ν)ℱ^2/ν,

as well as

    ν∫_0^t∇(s)^2 ds ≤_0^2+([t]+1)ℱ^2/ν,

where C_ denotes the  constant in the Poincaré inequality ^2≤
C_Ω∇^2, valid for all ∈ V.


  We focus on the proof of the a priori estimate (<ref>), the
  estimate (<ref>) being a direct consequence of the energy balance.
   By the Young inequality we deduce from the the energy inequality, 

   
    [                                        ∀ ξ, τ∈_+  s.t.   0 ≤ξ≤τ,;                                                                 ; (τ)^2+
           ν∫_ξ^τ∇(s)^2 ds≤(ξ)^2+
           1/ν∫_ξ^τ(s)^2_V' ds. ]

In particular, when 0 ≤τ - ξ≤ 1,  


    (τ)^2+
        ν∫_ξ^τ∇(s)^2 ds≤(ξ)^2+
      ℱ^2/ν.


From this point, we argue step by step. The case 0 ≤ t ≤ 1 is the first step, which
is straightforward. The second step is the heart of the proof.  The issue is that energy
may increase, without control, when the time increases.

We will show that even if this happens, we can still keep the control on it, thanks
to (<ref>). The last step is the concluding step, carried out by
induction on n writing t = τ +n, for τ∈ [0,1].



Step 1.  t ∈ [0,1]: take  ξ= 0,  t = τ∈ [0,1]. Then 
  
    (t) ^2 ≤_0 ^2 + ℱ^2/ν    ∀ t∈[0,1].

  

Step 2.  We show that the following implication holds true


    (t ) ≤(t+1 )  ⇒ {[   (t) ≤(1+C_/ν)ℱ^2/ν,; (t+1) ≤(2+C_/ν)ℱ^2/ν. ].


In the following we will set 

    𝒞^2:= 
      ( _0 ^2 + ℱ^2/ν ) +
    (2+C_/ν)ℱ^2/ν = _0^2+(3+C_/ν)ℱ^2/ν.



Sub-Step 2.1. By the energy inequality with ξ = t, and
τ = t+1, we have, by using the hypothesis on the L^2-norm at times t and t+1: 

    ν∫_t^t+1∇(s)^2 ds
           ≤(t+1)^2-
        (t)^2+ν∫_t^t+1∇(s)^2 ds
    
      
    ≤ℱ^2/ν.
 
Hence, by the Poincaré's inequality:


    ∫_t^t+1(s)^2 ds ≤    C_∫_t^t+1∇(s)^2 ds ≤C_ℱ^2/ν^2,


Sub-Step 2.2.  Let be given ϵ>0 and let ξ∈
[t,t+1] be such that

 
    (ξ)^2<inf_s∈[t,t+1](s)^2+ϵ≤(s)^2+ϵ    ∀ s∈ [t,t+1].


Let us write:
 
    (t)^2≤(t+1)^2   =(t+1)^2- (ξ)^2 +
        (ξ)^2
       =  (t+1)^2- (ξ)^2 + ∫_t^t+1(ξ)^2 ds,

being the integration with respect to the s variable.

To estimate the right-hand side we use the energy inequality with τ=t+1 to get


    (t+1)^2- (ξ)^2≤ℱ^2/ν.


Moreover, using the estimate (<ref>) we get,
 
    ∫_t^t+1(ξ)^2 ds ≤∫_t^t+1 (
        (s)^2+ϵ) ds ≤C_ℱ^2/ν^2 + ,

    therefore, letting → 0 yields, 
 
    (t)^2≤(1+C_/ν)ℱ^2/ν.

In addition, we get from the energy inequality

    (t+1 )^2≤(t) ^2 + ℱ^2/ν≤(2+C_/ν)ℱ^2/ν.


Step 3. Conclusion of the proof of (<ref>). Any
t≥0 can be decomposed as 

    t =n  + τ,   with n∈ and τ∈ [0,1[.


We argue by induction on n.  If n = 0, estimate (<ref>) has been proved
in Step 1. Assume that (<ref>) is satisfied for  t:=n+τ, for all n≤ N that
is, 

    (n + τ ) ≤𝒞     n=0,…,N,

where the constant 𝒞 is defined
by (<ref>).

If (N +1+ τ ) <  (N + τ ), then (<ref>)
holds at the time t=N+1+τ, by the inductive hypothesis.


If   (N + τ ) ≤ (N + 1+ τ ), then the
inequality (<ref>) is satisfied by  Step 2 for t=N+1+τ, ending the proof.  




Once we have proved that the uniform (independent of n∈) L^2-estimate is satisfied
by the Galerkin approximate functions, it is rather classical to prove that we can extract
a sub-sequence that converges weakly* in L^∞(0,T;L^2(Ω)) (for all positive
T) to a weak solution to the NSE, which inherits the same bound. We refer to the
references already mentioned for this point.



§ REYNOLDS DECOMPOSITION AND TIME-AVERAGING 



We sketch the standard routine, concerning time-averaging, when used in turbulence
modeling practice.  In particular, we recall the Reynolds decomposition and the Reynolds
rules. Then, we give a few technical properties of the time-averaging operator M_t, defined by


    M_t (ψ) := 1/t∫_0^t ψ (s)  ds,

for a given fixed time t > 0.  We need to apply it not only to real functions of a real
variable, but also to Banach valued functions, hence we need to deal with the Bochner integral.

Before all, we start with the following corollary of
Bochner theorem (see Yosida <cit.>).


   
  Assume that, for some t>0 we have ψ∈ L^p([0,t]; X) (namely ψ is a Bochner
  p-summable function over [0,t], with values in the Banach space X). Then, it holds 
  
    M_t(ψ)_X≤1/t^1/pψ_L^p([0,t];X).


Estimate (<ref>) is the building block to give a sense to the long-time
average as 

    ψ:= lim_t →∞ M_t (ψ),

whenever the limit exists.

It is worth noting at this stage that the mapping μ behind the long-time average,
defined on the Borel sets of _+ by
 
    A↦lim_t→+∞1/tλ(A∩[0,t])
       :=lim_t→+∞ M_t ( 1_A ) = μ(A),

 where λ the Lebesgue measure, 
 is not –strictly speaking– a probability measure since it is not
 σ-additive[The mapping μ satisfies μ(A∪ B) = μ(A) + μ (B)
   for A ∩ B = ∅ but, on the other hand, we have ∑_n=0^∞μ([n,
   n+1[) = 0 ≠ 1 = μ ( ⋃_n=0^∞ [n, n+1[ ). ]. Therefore, the quantity
 ψ is not rigorously a statistic, even if practitioners could be tempted to
 write it (in a suggestive and evocative meaningful way) as follows:
 
    ψ() = ∫__+ψ(s, )  d μ(s).




 §.§ General setup of turbulence modeling


We recall that M_t is a linear filtering operator which commutes with differentiation
with respect to the space variables (the so called Reynolds rules). In particular, one has
the following result (its proof is straightforward), which is essential for our modeling
process.

  Let be given ψ∈ L^1([0,T], W^1,p(Ω)), then 
  
    D M_t(ψ)=   M_t(Dψ)    ∀ t >0,

  for any first order differential operator D acting on the space variables
  ∈Ω. 


By denoting the long-time average of any field ψ by ψ as
in (<ref>), we consider the fluctuations ψ' around the mean
value, given by the Reynolds decomposition

    ψ :=ψ+ψ'.
 
Observe that long-time averaging has many convenient formal mathematical properties, recalled in the following.

 
  
  The following formal properties holds true, provided the long-time averages do exist.
  
  
  * The bar operator preserves the no-slip boundary condition. In other words,
    if ψ_|_Γ= 0, then ψ_|_Γ =0;

  
  * Fluctuation are in the kernel of the bar operator, that is ψ'=0;

  
  * The bar operator is idempotent, that is ψ=ψ, which also
    yields ψ φ=ψφ.
  


Accordingly, the velocity components can be decomposed in the Reynolds
decomposition as follows:

    (t, ) = () + '(t, ).

Let us determine (at least formally) the equation  satisfied by . To do so, we use the
above Reynolds rules to expand the nonlinear quadratic term into

    ⊗ =
    ⊗+'⊗',
 
which follows by observing that ' ⊗ = ⊗' = 0.  

The above rules allow us to prove the following result showing a certain “orthogonality”
between averages and fluctuations.

  
  Let be given a linear space X⊆ L^2(Ω) with a scalar product
( . , . ). Let in addition be given a function ψ: ^+→ X such that
ψ is well defined. Then it follows that
  
    (ψ',ψ')=(ψ,ψ)-(ψ,ψ),

provided all averages are well defined.


  The proof follows by observing that ψ'=ψ-ψ, hence 
  
    (ψ',ψ')=
          (ψ-ψ,ψ-ψ)=(ψ,ψ)-
        2(ψ,ψ)+(ψ,ψ),

and by the Reynolds rules (ψ,ψ)=(ψ,ψ) and
(ψ,ψ)=(ψ,ψ),
from which is follows the thesis.

In particular, we will use it for the V scalar product showing that 

    ∇'^2=∇^2-∇^2,

for  : ^+→ V, such that the long-time average exists. 


  Observe that, for weak solutions of the NSE , the average
  M_t(∇^2) is bounded uniformly, by the result of
  Theorem <ref> and –up to sub-sequences– some limit can be
  identified. Moreover, by using an argument similar to Lemma <ref> it follows
  that
  
    M_t(∇)^2≤ M_t(∇^2),

  which show that (up to sub-sequences) also the second term from the right-hand side can
  be properly defined. Consequently, also the average of the squared V-norm of the
  fluctuations from the left-hand side is well defined by difference.

Long-time averaging applied to the Navier-Stokes equations (in a strong
formulation) gives the following “equilibrium problem” for the long-time average
(),


    {
          -νΔ+∇·(⊗)+∇p   =    in Ω,
          
    ∇·   =0    in Ω, 
          
       =0    on Γ,
        .


which we will treat in the next section to make appear the closure problem. 


 §.§ Reynolds stress and Reynolds tensor

The first equation of system (<ref>) can be rewritten also as follows
(by using the decomposition into averages and fluctuations)


    -νΔ+∇·(⊗)+∇p=
      -∇·('⊗')+,

called the Reynolds equations.  Beside convergence issues, a relevant point is to
characterize the average of product of fluctuations from the right-hand side, which is the
divergence of the so called Reynolds stress tensor, defined as follows

    ='⊗'.
 
The Boussinesq hypothesis, formalized in <cit.> (see also <cit.>, for a comprehensive and modern presentation) corresponds then to a closure 
hypothesis with the following linear constitutive equation:

    =-ν_t∇+∇^T/2 + 2/3 k  Id,
 
where ν_t≥0 is a scalar coefficient, called turbulent viscosity
or eddy-viscosity (sometimes called “effective viscosity”), and

    k = 1/2| ' |^2,

is the turbulent kinetic energy,
see <cit.>. Formula (<ref>) is a linear relation between
stress and strain tensors, and shares common formal points with the linear constitutive
equation valid for Newtonian fluids.  In particular, this assumptions motivates the fact
that  must be dissipative[The sign adopted in (<ref>)
  is a convention consistent with our mathematical approach. However, according to the
  analogy of the Reynolds stress with viscous forces, it is also common to define it as :=- '⊗', which does not change anything.] on the mean
flow. Some recent results in the numerical verification of the hypothesis can be found in
the special issue <cit.> dedicated to Boussinesq.  Here, we show that, beside the
validity of the modeling assumption (<ref>), the Reynolds stress tensor
 is dissipative, under minimal assumptions on the regularity of the data of the
problem.



 §.§ Time-averaging of uniformly-local  fields
 

We list in this section some technical properties of the operator M_t acting on
uniform-local fields, and the corresponding global weak solutions to the NSE.  The first
result is the following
 
   
  Let 1 < p < ∞ and let be given f ∈ L^p_uloc(_+;X). Then
  
    ∀  t ≥ 1,      M_t (f) _X ≤ 2 f_L^p_uloc(_+;X).
 

 
  Applying (<ref>) and some straightforward inequalities yields
  
  
    M_t(f)_X≤1/t∫_0^t f _X ds
        ≤1/t∫_0^[t]+1f_X ds 
        ≤1/t∑_k=0^[t]∫_k^k+1 f _X ds.

  
  Therefore by the Hölder inequality we get:

  
    M_t(f)_X   ≤1/t∑_k=0^[t](∫_k^k+1 f _X^p ds)^1/p(∫_k^k+11 ds)^1/p'
       ≤[t]+1/t f _L^p_uloc(X)^p ≤2  f
          _L^p_uloc(_+;X)^p,


the last inequality being satisfied since [x]+1/x≤ 2 is valid for all
x≥1.


In the next section, we will focus on the case p=2 and X = V'. We will need the
following result, which is a consequence of Lemma <ref>.


   
  Let be given _0∈ H, f∈ L^2_uloc(_+;V'), and let  be a
  global weak solution to the NSE corresponding to the above data. Then, we have, ∀ t≥1, 

 
    M_t(∇^2) ≤_0^2/ν
          t+2ℱ^2/ν^2, 
            M_t( _V' ^2 ) ≤ 2 ℱ^2,
 

  where F =  f_L^2_uloc(_+;V').



  It suffices to divide estimate (<ref>) by ν t.  Therefore, it follows that


    1/t∫_0^t∇(s)^2 ds≤_0^2/ν t+2ℱ^2/ν^2.

Estimate (<ref>)  is straightforward. 

In particular the family { M_t () }_t ∈_+ is bounded in V. Therefore, we
can as of now state the following, which will be recalled in a more precise form in the
next sections.
 
  There exists  such that -up to a sub-sequence- M_t () → as t→∞,
  in appropriate topologies.
 
The following result is a direct consequence of (<ref>)
and (<ref>) combined with Cauchy-Schwarz inequality.

 
  The family { M_t( <  ,  > ) }_t >0 is bounded uniformly in t, and one has
  
    | M_t( <  ,  > ) | ≤√(2)ℱ ( _0^2/ν
            t+2ℱ^2/ν^2 )^1/2.

 

We finish this section with a last technical result, that we will need to prove 
Item <ref>)  of Theorem <ref>.


   
  Let 1 < p < ∞ and let be given f ∈ L^p_uloc(_+;X), which satisfies in
  addition
  
    ∃ f∈ X,   such that  lim_t→+∞∫_t^t+1 f(s)-f_X^p ds=0.
 
  Then, we have
  
    lim_t→+∞1/t∫_0^t f (s) -  f _X^p ds=0.
 
Moreover, M_t(f) weakly converges to f in X when t →+∞. In
particular, we have f = f.


 
  By the hypothesis (<ref>), we have that 

    ∀ ε>0  ∃  M∈:  ∫_t^t+1f (s) -
      f_X^p ds<ε/2    ∀ t>M.

Hence, for  t ≥ M, then 

    1/t∫_0^tf(s) -
          f_X^p ds=1/t∫_0^Mf (s) -
          f_X^p ds+1/t∫_M^tf(s) -
          f_X^p ds
        
    ≤M/t (  f _L^p_uloc(X)^p+ f _X^p )
        +[t ]+1-M/tε/2.

It follows that one can choose M large enough such that


    1/t∫_0^tf (s)- f _X^p ds<ε    ∀  t>M,


hence, being this valid for arbitrary ε>0, it
follows (<ref>). 

It remains to prove the weak convergence of M_t(f) to f∈ X when
t→+∞. To this end, let be given φ∈ X'. Then, we have


    < φ, M_t(f) > - < φ, f > = 1/t∫_0^t
      <φ, f(s)- f >  ds,


which leads to 


    | < φ, M_t(f) > - < φ, f > | ≤1/t∫_0^t  φ_X'  f(s)- f_X  ds,

and by Hölder inequality, 


    | < φ, M_t(f) > - < φ, f > | 
      ≤φ_X'(  1/t∫_0^t    f(s)- f_X^p  ds  )^1/p,


yielding, by (<ref>), to lim_t →+∞ < φ, M_t(f)> =<φ, f >, hence concluding the proof.
 
The following corollary definitively concludes Item <ref>) of
Theorem <ref>.
 
  Let be given _0∈ H, f∈ L^2_uloc(_+;V') that
  satisfies (<ref>), and let  be a global weak solution to the NSE
  corresponding to the above data. Moreover let  be such that lim_t →∞ M_t  = in V (eventually up to a sub-sequence), then


    lim_t →∞ M_t ( <  ,  > ) =   < ,  >.

 
 
 Let us write the following decomposition:
  
    1/t∫_0^t<,> ds=
        1/t∫_0^t<-
        ,> ds+1/t∫_0^t< ,> ds.

  On one hand since ∈ V is independent of t, we obviously have
  
    1/t∫_0^t<,> ds→  <,>.
 
  On the other hand, we have also 
  
    |1/t∫_0^t<- ,> ds
          |   ≤1/t∫_0^t- _V'∇ ds
          
       ≤(1/t∫_0^t-
            _V'^2 ds)^1/2(1/t∫_0^t∇^2 ds)^1/2.
 
  Combining (<ref>) with (<ref>) shows that the
  right-hand side in (<ref>) vanishes as t →∞.  

  It is important to observe that 
  
    M_t (<,> )→ <,>,

  is –in some sense– an assumption on the (long-time) behavior of the “covariance”
  between the external force and the solution itself.  Cf. Layton <cit.> for a
  related result in the case of ensemble averages.

  The control of the (average/expectation of) kinetic energy in terms of the energy input
  is one of the remarkable features of classes of statistical solutions, making the
  stochastic Navier-Stokes equations very appealing in this context. See the review, with
  applications to the determination of the Lilly constant, in Ref. <cit.>. See
  also <cit.>.




§ PROOF OF THEOREM <REF>



In all this section we have as before _0 ∈ H, f∈ L^2_uloc(_+;V'),
and  is a global weak solution to the NSE (<ref>) corresponding to the above
data.  We split the proof of Theorem <ref> into two steps. We first apply
the operator M_t to the NSE, then we extract sub-sequences and take the limit in the
equations. In the second step we make the identification with the Reynolds stress 
and show that it is dissipative in average, at least when f satisfies in
addition (<ref>).



 §.§ Extracting sub-sequences
 
We set:

    _t():=M_t()().

Applying the operator M_t on the NSE we see that for almost all t≥0 and for all
∈ V, the field _t is a weak solution of the following steady Stokes problem
(where t>0 is simply a parameter)

    ν∫_Ω∇_t:∇ d+∫_ΩM_t((·∇) )· d   = < M_t(),>
        
       + ∫_Ω_0-(t)/t· d.


The full justification of the equality (<ref>) starting from the
definition of global weak solutions can be obtained by following a very well-known path
used for instance to show with a lemma by Hopf that Leray-Hopf weak solutions can be
re-defined on a set of zero Lebesgue measure in [0,t] in such a way that (s)∈ H
for all s∈[0,t], see for instance Galdi <cit.>. In fact,
by following ideas developed among the others by Prodi <cit.>, one can take
χ_[a,b] the characteristic function of an interval [a,b]⊂, and use as
test function its regularization multiplied by ∈ V. Passing to the limit as the
regularization parameter vanishes one gets (<ref>).

The process of extracting sub-sequences, which is the core of the main result, is reported
in the following proposition.


  
  Let be given a global solution  to the NSE, corresponding to the data _0∈
  H and ∈ L^2_uloc(_+;V'). Then, there exist
  
  
  * a sequence {t_n}_n∈ that goes to +∞
    when n goes to +∞;
  
  * a vector field ∈ V'; 
 
  * a vector field ∈ V;
 
  * a vector field ∈ L^3/2(Ω)^3;  

such that such that it holds when n →∞:
   
    M_t_n()⇀    in  V',
           
    
              M_t_n()⇀    in  V,
           
       M_t_n ((·∇)  )⇀    in  L^3/2(Ω)^3⊂ V',

and for all ∈ V

    ν∫_Ω∇:∇ d
      +∫_Ω· d= <,>.

Moreover, by defining

    :=-(·∇) ∈ L^3/2(Ω)^3,
 
we can also rewrite (<ref>) as follows

    ν∫_Ω∇:∇ d+
      ∫_Ω(·∇) · d
      +∫_Ω· d= < ,>;


  * writing ' = -, ' =  -, we also have


    M_t_n ( < ,  > ) →   <,  > +  < ', ' > .

      
  


  As f∈ L^2(_+;V'), we deduce from Lemma <ref> that
  {M_t( f)}_t>0 is bounded in V'.  Hence, we can use weak pre-compactness of
  bounded sets in the Hilbert space V' to infer the existence of t_n and
  ∈ V' such that M_t_n()⇀ in V'.

  Next, estimate (<ref>) from Lemma <ref>,
  combined with estimate (<ref>) from Lemma <ref>, leads to the bound
  
    ∃ c>0:    ∇ M_t() = M_t(∇)≤ c     ∀ t>0,

  proving (up to a the extraction of a further sub-sequence from {t_n}, which we call
  with the same name) that M_t_n()⇀ in V'.

  Then, we observe that, if ∈ L^∞(0,T;H)∩ L^2(0,T;V) by classical
  interpolation


    (  ·∇) ∈ L^r(0,T;L^s(Ω))    with
        2/r+3/s=4,    r∈[1,2].

In particular, we get

    (
      ·∇) _L^3/2(Ω)≤_L^6∇_L^2≤ C_S∇^2,

where C_S is the constant of the Sobolev embedding H^1_0(Ω)↪
L^6(Ω). Hence, by using the bounds (<ref>)-(<ref>) on
the weak solution  we obtain that


    ∃ c:    M_t ((·∇) ) )_L^3/2(Ω)≤
      c,    ∀ t>0,

proving that, up to a further sub-sequence relabelled again as {t_n},

    M_t_n ((·∇)  )⇀    in  L^3/2(Ω)^3,

for some vector field ∈ L^3/2(Ω)^3.

Next, we use (<ref>) which shows that

    ∫_Ω_0-(t)/t· d→0    as t→+∞.

Then, writing the weak formulation and by using the results of weak convergence previously
proved, we get (<ref>). Then, the identity (<ref>) comes simply
from the definition (<ref>) of .

It remains to prove the last item. We know from (<ref>) that the sequence { M_t_n (< ,  > ) } _n ∈ is bounded in . By extracting again a
sub-sequence (still denoted by t n), we can get a convergent sequence still
denoted (after relabelling) by { M_t_n (< ,  > ) } _n ∈, and let
< ,  > be its limit. Let us write the decomposition


    [                                                               M_t_n (< ,  >);                                                                             ; 1cm 
        =<,  > + M_t_n (< ',  >) + M_t_n (< , ' >)+
        M_t_n ( < ', ' > ). ]


As M_t_n (< ',  >) =   < M_t_n ( '),  >, we deduce from the results
above that M_t_n (< ',  >) → 0 as n →∞. Similarly, we also have
M_t_n (< , ' >) → 0. Hence, we deduce
from (<ref>) that { M_t_n ( < ', ' > ) }_n ∈ is
convergent, and if we denote by < ', ' > its limit, the following
natural decomposition holds true:

    < ,  > = <,  >+< ', ' >,

concluding the proof.




 §.§ Reynolds stress, energy balance and dissipation



In the first step we have identified a limit (,) for the
time-averages of both velocity and external force (,). We need now to recast this
in the setting of the Reynolds equations, in order to address the proof of the Boussinesq
assumption.



  Beside the results in Proposition <ref>, in order to complete the proof of
  Theorem <ref>, we have to prove the following facts:

  
  
  *  the proper identification of the limits with the Reynolds stress
    ;

  
  *  the energy balance for ;

  
  *  to prove that (<ref>)  holds, namely
  
    =  ν∇'^2≤∫_Ω(∇·)· d + < ', ' > .
  


We proceed in the same order.

Item <ref>.  Since ∈ L^2(0,T;V)⊂ L^2(0,T;L^6(Ω)^3), it
  follows that ⊗∈ L^1(0,T;L^3(Ω)). Hence, the same argument as in the
  previous subsection shows that (possibly up to the extraction of a further sub-sequence)
  there exists a second order tensor θ∈ L^3(Ω)^9 such that
  
    M_t_n(⊗)⇀θ    in L^3(Ω)^9.

  Let us set
  
    :=θ-⊗.

  Since the operator M_t commutes with the divergence operator, the
  equation (<ref>) becomes
  
    ν∫_Ω∇_t:∇ d- ∫_ΩM_t(⊗):  d   = < M_t(),>
          
       + ∫_Ω_0-(t)/t· d.

  Then, by taking the limit along the sequence t_n→+∞
  in (<ref>), we get[
    According to the formal decomposition (<ref>), this
      suggests that M_t_n(' ⊗') → 0, provided that one is able to give
      a rigorous sense and sufficiently strong bounds on '⊗', for the weak
      solution .


] the equality
  
    =∇·.

  

  Item <ref>.  We use ∈ V in (<ref>) as test function
  and we obtain the equality
  
    ν∇^2+∫_Ω
        (∇·)· d= <,>.
  
  We observe that due to the absence
  of the time-variable the following identity concerning the integral over the space
  variables is valid
  
    ∫_Ω (·∇) · d=
        ∫_Ω (·∇) ||^2/2 d=0    ∀ ∈ V.


  This is one of the main technical facts which are typical of the mathematical analysis
  of the steady Navier-Stokes equations and which allow to give precise results for the
  averaged Reynolds equations. On the other hand, we recall that if (t,) is a
  non-steady (Leray-Hopf) weak solution, then the space-time integral 
  
    ∫_0^T∫_Ω
        (·∇) · d dt,

  is not well defined and consequently the above integral vanishes only formally.


  Item <ref>. From now, we assume that the
  assumption (<ref>) in the statement of
  Theorem <ref> holds true. We integrate the energy
  inequality (<ref>) between 0 and t_n and we divide the result by t_n>0,
  which leads to
  
    (t)^2/2t_n+1/t_n∫_0^t_n∇(s)^2 ds≤_0^2/2t_n+1/t_n∫_0^t_n<,> ds.

  Recall that by Lemma <ref>
  
    (t)^2/2t→0  and  _0^2/2t→0  as t→+∞.

   Therefore, we take the limit in (<ref>) and we use (<ref>), which yields
 
    ν^2≤<  ,  > =    <,  > +  < ', ' >.

By (<ref>)  we then have

    ν^2≤ν∇^2+∫_Ω
        (∇·)· d +   < ', ' >,

  which yields (<ref>) by (<ref>), concluding the
  proof.
 
 



§ ON ENSEMBLE AVERAGES



In this section we show how to use the results of Theorem <ref> to give
new insight to the analysis of ensemble averages of solutions. In this case we study
suitable averages of the long-time behavior and not the long-time behavior of statistics,
as in Layton <cit.>.

Since we first take long-time limits and then we average the Reynolds
equations, the initial datum is not so relevant. In fact due to the fact that it holds

    _0^2/t→0    as t→+∞,

then the mean  is not affected by the initial datum.

As claimed in the introduction, we consider now the problem of having several external
forces, say a whole family {^k}_k∈⊂ V', all independent of time. We
can think as different experiments with slightly different external forces, whose
difference can be due to errors in measurement or in the uncertainty intrinsic in any
measurement method. In particular, one can consider for a given force  and that
{^k} will represent small oscillations around it, hence we can freely assume that
we have an uniform bound


    ∃ C>0:    ^k_V'≤ C    ∀ k∈.
 

Having in mind this physical setting, we denote by ^k∈ V the long-time
average of the solution corresponding to the external force ^k∈ V' and, as
explained before (without loss of generality) to the initial datum _0=0. The
vector ^k∈ V satisfies for all ∈ V the following equivalent
equalities for all k∈


    ν∫_Ω∇^k:∇ d+
      ∫_Ω^k· d= < ^k,>,
      
       ν∫_Ω∇^k:∇ d+
      ∫_Ω(^k·∇) ^k· d
      +∫_Ω^k· d= < ^k,>,

for appropriate ^k,^k∈ L^3(Ω)^3/2. Since both V and V' are Hilbert
spaces, by using (<ref>) it follows that there exists ∈ V'
and a sub-sequence (still denoted by {^k}) such that


    ^k⇀     in V'.


Our intention is to characterize, if possible, the limit of
{^k}_k∈. If the forces are fluctuations around a mean value, then the
field ^k will remain bounded in V, but possibly without converging to some
limit. From an heuristic point of view one can expect that averaging the sequence of
velocities (which corresponds to averaging the result over different realizations) one can
identify a proper limit, which retains the “average” effect of the flow. 

Again, it comes into the system the main idea at the basis of Large Scale methods: The
average behavior of solutions seems the only quantity which can be measured or
simulated. 

It is well-known that one of the most used summability technique is that of
Cesàro and consists in taking the mean values, hence we focus on the arithmetic mean of
time-averaged velocities

    𝐒^n:=1/n∑_k=1^n ^k.

It is a basic calculus result that if a real sequence {x_j}_j∈ converges to
x∈, then also its Cesàro mean S_n = 1/n∑_j=1^n x_j will converge to
the same value x. On the other hand, the converse is false; sufficient conditions on the
sequence {x_j}_j∈ implying that if the Cesàro mean converges, then the
original sequence converges, are known in literature as Tauberian theorems. This is a
classical topic in the study of divergent sequences/series.  In the case of X-valued
sequences {^k}_k∈ (the space X being an infinite dimensional Banach space)
one has again that if a sequence converges strongly or weakly, then its Cesàro mean will
converge to the same value, strongly or weakly in X, respectively.

The fact that averaging generally improves the properties of a sequence, is reflected also
in the setting of Banach spaces even if with additional features coming into the theory.
Two main results we will consider are two theorems known as Banach-Saks and Banach-Mazur.

Banach and Saks originally in 1930 formulated the result in
L^p(0,1), but it is valid in more general Banach spaces.


  Let be given a bounded sequence {x_j}_j∈ in a reflexive Banach space
  X. Then, there exists a sub-sequence {x_j_k}_k∈ such that the sequence
  {S_m}_m∈ defined by
  
    S_m := 1/m∑_k=1^m x_j_k,

converges strongly in X.

The reader can observe that in some cases it is not needed to extract a
sub-sequence (think of any orthonormal set in an Hilbert space, which
is weakly converging to zero, and the Cesàro averages converge to
zero strongly), but in general one cannot infer that the averages of
the full sequence converge strongly. One sufficient condition is that
of uniform weak convergence. We recall that {x_j}⊂ X
uniformly weakly converges to zero if for any ϵ>0
there exists j∈, such that for all ϕ∈ X', with
ϕ_X'≤1, it holds true that


    #{ j∈ N:  |ϕ(x_j)|≥ε}≤ j.


See also Brezis <cit.>.


Another way of improving the weak convergence to the strong one is by the by the
convex-combination theorem (cf. Yosida <cit.>).


  Let (X, . _X) be a Banach space and let {x_j}⊂ X be a sequence
  such that x_j⇀ x as j→+∞.

  Then, one can find for each n∈, real coefficients {α^n_j}, for
  j=1,…,n such that
  
    α^n_j≥0    and     ∑_j=1^nα^n_j=1,

  such that 
  
    ∑_j=1^nα^n_j x_j→ x    in X,  as  
        n→+∞,

  that is we can find a “convex combination” of {x_j}, which strongly converges to
  x∈ X.











One basic point will be that of considering averages of the external forces, which we will
denote by ^n and considering the same averages of the solution of the
Reynolds equations ^n.  They are both bounded and hence, weakly
converging (up to a sub-sequence) to ∈ V' and ∈ V,
respectively. Then, in order to prove that the dissipativity is preserved one has to
handle the following limit of the products

    lim_n→+∞  <^n,^n>,

which cannot be characterized, unless (at least) one of the two terms
converges strongly. This is why we have to use special means instead
of the simple Cesàro averages

The first result of this section is then the following:


  
  Let be given {^k}_k∈ uniformly bounded in V'. Then one can find either a
  Banach-Saks sub-sequence or a convex combination of {^k}_k∈, which are
  converging weakly to some ∈ V, which satisfies a Reynolds
  system (<ref>), with an additional dissipative term.


  We define ^n and ^n to be either 
  
    ^n:=1/n∑_k=1^n^j_k    and    ^n:=1/n∑_k=1^n^j_k,

  or alternatively
  
    ^n:=∑_j=1^nα^n_j^j    and    ^n:=∑_j=1^nα^n_j^j,

  where the sub-sequence {j_k}_k∈ or the coefficients
  {α^n_j}_j, n∈ are chosen accordingly to the Banach-Saks of Banach-Mazur
  theorems in such a way that in both cases

    ^n→    in V'.



We define, accordingly to the same rules ^n, and we observe
that, by linearity, we have ∀ n∈
  
    ν∫_Ω∇^n:∇ d+
        ∫_Ω^n· d= < 
        ^n,>    ∀ ∈ V.


  Then, we can define
  ^n:=^n-(^n·∇) ^n, to
  rewrite (<ref>) also as follows

  
    ν∫_Ω∇^n:∇ d+
        ∫_Ω(^n·∇) ^n· d+∫_Ω^n· d
        =  <^n,>.

  By the uniform bound on ^k_V' and by results of
  Section <ref> on the Reynolds equations it
  follows that there exists C such that ^k_V≤ C,
  hence 
  
    ^n_V  ≤ C    ∀ n∈,

  and  we can suppose that (up to sub-sequences) we have weak
  convergence of the convex combinations
  
    ^n   ⇀     in V,
          
    ^n   ⇀     in L^3/2(Ω)^3,
          
    ^n   ⇀     in L^3/2(Ω)^3,

  Hence, passing to the limit in (<ref>), we obtain
  
    ν∫_Ω∇:∇ d+
        ∫_Ω· d= < 
        ,>    ∀ ∈ V.


  By the same reasoning used before we have, for
  :=-(·∇),

  
    ν∫_Ω∇:∇ d+
        ∫_Ω(·∇) · d+
        ∫_Ω· d= <  
        ,>.


  Then, if we take = in (<ref>) we obtain

  
    ν∇^2
        +∫_Ω· d= <,>.


  On the other hand, if we take =^n in (<ref>) and by the
  result of the previous section, we have

  
    ν∇^n^2≤ <^n,^n>,

  hence passing to the limit, by using the strong convergence of
  ^n in V' and the weak convergence of ^n
  in V we have
  
    ν∇^2≤lim inf_n→+∞ν∇^n^2≤ 
        <,>.

  If we compare with (<ref>) we have finally the dissipativity

  
    1/|Ω|∫_Ω(∇·)· d
        =1/|Ω|∫_Ω· d≥0,

  that is a sort of ensemble/long-time Boussinesq
  hypothesis, cf. with the results from Ref. <cit.>.


In the previous theorem, we have a result which does not concern directly with the
ensemble averages, but a selection of special coefficients is required. This is not
completely satisfactory from the point of view of the numerical computations, where the
full arithmetic mean should be considered. The main result can be obtained at the price of
a slight refinement on the hypotheses on the external forces

To this end we recall a lemma, which is a sort of Rellich theorem in
negative spaces (see also Galdi <cit.> and
Feireisl <cit.>).



Let Ω⊂^n be bounded and let be given 1<p<n. Let {f_k}_k∈ be a
sequence uniformly bounded in L^q(Ω) with q>(p^*)', where p^*=n p/n-p is
the exponent in the Sobolev embedding W^1,p_0(Ω)↪
L^p^*(Ω). Then, there exists a sub-sequence {f_k_m}_m∈ and f∈
L^q(Ω) such that
  
    f_k_m⇀ f    in L^q(Ω),
         
        f_k_m→ f    in W^-1,p'(Ω),

  or, in other words, the embedding L^q(Ω)↪↪
  W^-1,p'(Ω) is compact.


We present the proof for the reader's convenience.

 
Since by hypothesis L^q(Ω) is reflexive, by the
Banach-Alaouglu-Bourbaki theorem we can find a sub-sequence  f_k_m
such that 

    f_k_m⇀ f    in L^q(Ω),

and by considering the sequence { f_k_m- f}_m∈ we can suppose that f=0.  We
then observe that by the Sobolev embedding we have the continuous embeddings

    W^1,p_0(Ω)↪  L^p^*(Ω)∼
      L^(p^*)'(Ω)↪  (W^1,p_0(Ω))'≃
      W^-1,p'(Ω),

where L^p^*(Ω)∼ L^(p^*)'(Ω) is the duality identification, while the
second one (W^1,p_0(Ω))'≃ W^-1,p'(Ω) is the Lax isomorphism. This
shows L^(p^*)'(Ω)⊆ W^-1,p'(Ω).

Next, let be given a sequence {f_k_m}⊂ W^-1,p'(Ω), then
by reflexivity (since 1<p<∞) there exists {ϕ_k_m}⊂
W^1,p_0(Ω) such that

    f_k_m_W^-1,p'(Ω)=f_k_m(ϕ_k_m)= <
        f_k_m,ϕ_k_m>,
 with
  ϕ_k_m_W^1,p_0(Ω)=∇ϕ_k_m_L^p(Ω)=1.

Hence, by using the classical Rellich theorem, we can find a
sub-sequence {ϕ_k_j}_j∈ such that

    ϕ_k_j→ϕ    in L^r(Ω)    ∀ r<p^*.

In particular, we fix r=q' (observe that q>(p^*)' implies
q'<p^*) and we have 

    f_k_m_W^-1,p'(Ω)= < f_k_m,ϕ_k_m-ϕ>+< f_k_m,ϕ>.

The last term converges to zero, by the definition of weak convergence
f_k_m⇀0, while the first one satisfies

    |<    f_k_m,ϕ_k_m-ϕ>|≤f_k_m_W^-1,p'ϕ_k_m-ϕ_W^1,p_0,

and since f_k_m_W^-1,p' is uniformly bounded and
ϕ_k_m-ϕ_W^1,p_0 goes to zero, then also this one vanishes as
j→+∞.




  The proof of this theorem can be obtained by following the same ideas of the
  Proposition <ref>.  In fact, the main improvement is that the weak
  convergence ^j⇀ in L^q(Ω) implies (without
  extracting sub-sequences) that
  
    ^k→    in V'.

  This follows since from any sub-sequence we can find a further
sub-sequence which is converging strongly, by
Lemma <ref>. Then, by the weak convergence of the original
sequence, the limit is always the same and this implies that the whole
sequence {^k} strongly converges to its weak limit.

Hence, we have

    1/n∑_k=1^n   ^k→    in V',

and then, since ^n⇀ in V, we can infer that

    < ^n,^n > =  <1/n∑_k=1^n   ^k,
      1/n∑_k=1^n^k> →   < , >,

and the rest follows as in Proposition <ref>.







§ ACKNOWLEDGMENTS

The research that led to the present paper was partially supported by a grant of the group
GNAMPA of INdAM the project the University of Pisa within the grant
PRA_2018_52 UNIPI Energy and regularity: New techniques for classical PDE
problems.

#10=#10=0 0
  by1ptto00
  .2em∘#1"17 #1
  ' #10=#11.5ex`0 '
10

2AG91 
C. Amrouche, and V. Girault. 
On the existence and Regularity of the solutions of Stokes Problem in arbitrary dimension. 
Proc. Japan Acad. Ser. A Math. Sci., 67(5):171–175, 1991. 






BF2010
L. C. Berselli and F. Flandoli.
On a stochastic approach to eddy viscosity models for turbulent
  flows.
In Advances in mathematical fluid mechanics, pages 55–81.
  Springer, Berlin, 2010.

BIL2006
L. C. Berselli, T. Iliescu, and W. J. Layton.
Mathematics of Large Eddy Simulation of turbulent flows.
Scientific Computation. Springer-Verlag, Berlin, 2006.


Bois2007
P.-A. Bois (Ed.).
Joseph Boussinesq, a scientist of bygone days and present times.
Comptes Rendus-Mécanique, Special vol. 335(9-10):475–678,  2007.


Bou18777
J. Boussinesq.
Theorie de l'ècoulement tourbillant.
Mém. prés par div. savants á la Acad. Sci., 23:46–50,
  1877.

Bre2011
H. Brezis.
Functional analysis, Sobolev spaces and partial differential
  equations.
Universitext. Springer, New York, 2011.

CL2014
T. Chacón Rebollo and R. Lewandowski.
Mathematical and numerical foundations of turbulence models and
  applications.
Modeling and Simulation in Science, Engineering and Technology.
  Birkhäuser/Springer, New York, 2014.

CF1988
P. Constantin and C. Foias.
Navier-Stokes equations.
Chicago Lectures in Mathematics. University of Chicago Press,
  Chicago, IL, 1988.



Fei2004
E. Feireisl.
Dynamics of viscous compressible fluids, volume 26 of Oxford Lecture Series in Mathematics and its Applications.
Oxford University Press, Oxford, 2004.

Fla2005
F. Flandoli.
An introduction to stochastic 3D fluids.
In SPDE in hydrodynamics: recent progress and prospects, volume
  1942 of Lecture Notes in Math., pages 51–148. Springer-Verlag, Berlin,
  2008.
Lectures from the C.I.M.E.-E.M.S. Summer School held in Cetraro,
  August 29-September 3, 2005.



Foi197273
C. Foias.
Statistical study of Navier-Stokes equations. I,
  II, Rend. Sem. Mat. Univ. Padova, 48:219–348, 1972; 
  ibid. 49:9–123, 1973.


FMRT2001a
C. Foias, O. Manley, R. Rosa, and R Temam.
Navier-Stokes equations and turbulence, volume 83 of Encyclopedia of Mathematics and its Applications.
Cambridge University Press, Cambridge, 2001.

UF95
U. Frisch.
Turbulence: The legacy of A. N. Kolmogorov
Cambridge University Press, Cambridge, 1995.

Gal2000a
G. P. Galdi.
An introduction to the Navier-Stokes initial-boundary value
  problem.
In Fundamental directions in mathematical fluid mechanics, Adv.
  Math. Fluid Mech., pages 1–70. Birkhäuser, Basel, 2000.

Gal2011
  G. P. Galdi.
  An introduction to the mathematical theory of the
    Navier-Stokes equations. Steady-state problems.
  Springer Monographs in Mathematics. Springer-Verlag, New York, 2011.
  
GR86
  V. Girault and P. A. Raviart.
  Finite element methods for Navier-Stokes equations.
  Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 1986. 
  
  
JL2016 
  N. Jiang and W. J. Layton.
  Algorithms and models for turbulence not at statistical equilibrium,
  Comput. Math. Appl.  71(11): 2352–2372, 2016.



LL2016
A. Labovsky and W. J. Layton.
Magnetohydrodynamic flows: Boussinesq conjecture.
J. Math. Anal. Appl., 434(2):1665–1675, 2016.

Lay2014
W. J. Layton.
The 1877 Boussinesq conjecture: Turbulent fluctuations are
  dissipative on the mean flow.
Technical Report TR-MATH 14-07, Pittsburgh Univ., 2014.


Ler1934 
J. Leray.
Sur le mouvement d'un liquide visqueux emplissant l'espace,
 Acta Math. 63(1), 193–248, 1934.


Lew2015
R. Lewandowski.
Long-time turbulence model deduced from the Navier-Stokes
  equations.
Chin. Ann. Math. Ser. B, 36(5):883–894, 2015.

JLL69 
  J.L. Lions.
   Quelques méthodes de résolution des problèmes aux limites non
    linéaires.
  Dunod, Paris, 1969. 
  
PL96 
  P.L. Lions.
   Mathematical topics in fluid mechanics. Vol. 1. Incompressible models. 
   Oxford Lecture Series in Mathematics and its Applications, 3. Oxford Science
    Publications. The Clarendon Press, Oxford University Press, New York, 1996.
  
  




 
 
  





  
LP25 
  L. Prandtl. 
  Über die ausgebildeten Turbulenz. 
  Z. Angew. Math. Mech., 5:136-139, 1925.  
  
  LP52 
  L. Prandtl. 
  Guide à travers la mécanique des fluides. 
  Dunod, Paris, 1952.  
  
Pro1959
  G. Prodi.
  Un teorema di unicità per le equazioni di Navier-Stokes,
  Ann. Mat. Pura Appl. (4) 48:173–182, 1959.

Rey1895
  O. Reynolds.
  On the dynamic theory of the incompressible viscous fluids and the
  determination of the criterion.
  Philos. Trans. Roy. Soc. London Ser. A, 186:123–164, 1895.
  

  
  
 
 
 

  2GS51 
 G.  Stokes. 
    On the Effect of the Internal Friction of Fluids on the Motion of Pendulums. 
      Transactions of the Cambridge Philosophical Society. 9:8–106 (1851)
  
LT06
  L. Tartar.
  An introduction to Navier-Stokes equation and oceanography. 
  Lecture Notes of the Unione Matematica Italiana,
  Springer-Verlag, Berlin; UMI, Bologna, 2006.
  
Tem1977b
  R. Temam.
  Navier-Stokes equations. Theory and numerical analysis.
  North-Holland Publishing Co., Amsterdam-New York-Oxford, 1977.
  
Yos1980
  K. Yosida.  
  Functional analysis.  
  Classics in Mathematics. Springer-Verlag, Berlin, 1995.
  Reprint of the sixth (1980) edition.
  



