RIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

SUB-GRID
TURBULENCE MODELLING
BY HOMOGENIZATION

Pierre PERRIER
Olivier PIRONNEAU

Rapport de Recherche N° 318

Août 1978
SUB-GRID TURBULENCE MODELLING BY HOMOGENIZATION

Pierre Perrier* and Olivier Pironneau**

ABSTRACT

If the intertial subrange of the turbulence of a detached turbulent flow is modelled by a stochastic
process satisfying the 1941-similarity hypothesis of Kolmogorov, the flow can be simulated by
solving numerically a Navier-Stokes equation with random coefficients. By approximating the
random coefficients by deterministic pseudo-periodic function we can use the theory of homoge-
nization to derive an equation for the mean flow. By this procedure we find that the mean flow sa-
tisfies a generalized Navier-Stokes equation with a random turbulent eddy viscosity tensor.

RESUME

Si le domaine inertiel de la turbulence d'un écoulement incompressible est modélisé par un proces-
sus stochastique satisfaisant aux hypothèses de similarité de Kolmogorov, alors on peut modéliser
l'écoulement complet par une équation de Navier-Stokes à coefficients pseudo-périodiques aléatoi-
res. On utilise alors la théorie de l'homogénéisation pour moyenner l'équation et on montre ainsi
que l'écoulement moyen satisfait à une équation de Navier-Stokes avec une viscosité tensorielle
turbulente.

* AMD-BA Industries, 75, Quai Carnot - 92210 ST CLOUD, France.

** IRIA-LABORIA, Domaine de Voluceau, Rocquencourt, 78150 LE CHESNAY, France.
# INTRODUCTION

Beside the analytical approaches, there are, at present, three main ways to compute turbulent flows (see for example MATHIEU (1977)). The first approach solves directly the full time dependent Navier-Stokes equations, either in the physical plane (REYNOLDS (1975)), or in the Fourier space (ORSZAG-PATTERSON (1971)) with velocity or vorticity as variables ; the limitation of this method in Reynolds number is obvious.

The second way consists of solving a time and space averaged of the Navier-Stokes equations ; a closure hypothesis for the Reynolds stresses is necessary as a substitute to the lack of information on the relation between the mean flow and the moments of the fluctuating velocities ; an eddy viscosity model is chosen in general (SMAGORINSKY (1963), DEARDORFF (1973), SCHUMANN (1975)).

The third way used in meteorology is to keep the full time dependent Navier Stokes equation on a coarse grid coupled with a subgrid scale turbulence simulation.

The present approach is of the last type. As for the Reynolds stresses closure problem, the modelling can be done in the physical plane or in the Fourier space depending on the type of equilibrium chosen ; the first will be described in terms of fluctuating velocities, momentum and time history of the momentum ; the second stresses spectral repartition of energy if the turbulence is isentropic. The subgrid scale modelling approach is unlikely to apply when many different scales are present in the flow and the turbulence is not in equilibrium at the subgrid size chosen.

In the present paper we assume that it is always possible to find at each point of the flow an eddy size such that the energy spectrum for larger wave numbers is of an equilibrium type, as measured experimentally (COMPTE-BELLOT and COURSIN (1970) and assumed theoretically (KOLMOGOROV (1941)) and such that at wave numbers two or more order of magnitude greater ther is no major dissipation. With the two previous assumptions it is possible to consider the turbulent subgrid flow as a well-defined random process which, when coupled with the main flow, perturbs it in a way that we shall evaluate by the theory of homogenization. This way the terms coupling the mean flow with the subgrid flow can introduce stochastic perturbations in the main flow, in addition to the creation of Reynolds
stresses and dissipation analogous to eddy viscosity effects.

If we analyse the turbulent flow past any complex body we can distinguish four main types. One corresponds to thin layers of turbulent equilibrium shear flows for which the greatest characteristic turbulent length is an order of magnitude smaller than the length of the turbulent layer and for which any parabolic approach with averaged Navier-Stokes equation will remain the best computational approach. Another one corresponds to an unsteady turbulent flow where it is possible to assume that the energy convected by large unsteady structures is some order of magnitude greater than the small scale turbulent energy so that the large structures are purely deterministic and coupled with the small structures by an eddy viscosity effect giving one predominant integral length scale to the turbulent flow. In that case all the classical subgrid modelling is to be used.

A third part of the flow is such that an equilibrium state has been reached in the spectral distribution of energy of small scale structures, due for example to an equilibrium of the first type of any upstream shear layer, but the large structures are either unstable or coupled to an external elliptic-type flow. In such cases the inverse cascade of energy is of main influence on the large structures of turbulent flow and the present theory can afford new computational possibilities.

The last type of flows that can be recognized correspond to a mixture of random and deterministic phenomenons of completely different scales in the same region of the flow. These cases are untractable with present computers though some rough mixing of the three different approaches could perhaps give some approximate simulation.

So the purpose of the paper is to derive an equation for the mean velocity $u^o$ of a turbulent flow $u$, whose turbulent part $u'$, is mainly homogeneous and isotropic. The classical approach for such flows is to write the solution to the Navier-Stokes equation :

$$ (1) \quad \begin{cases} \frac{\partial u}{\partial t} - v\Delta u + u\nabla u = \nabla p + f, \quad \nabla \cdot u = 0 \quad \text{in } \Omega \times ]0,T[ \\ u(x,o) = u_o; \quad u|_{\Gamma} = u_l \end{cases} $$
in the form

$$
(2) \quad u = u^0 + u'
$$

where $u'$ has zero mean value so that an averaging of (1), denoted by $\bar{\cdot}$, yields

$$
(3) \quad \begin{cases}
\displaystyle \frac{\partial u^0}{\partial t} - v \Delta u^0 + \overline{u^0 \nabla u'} + \overline{u' \nabla u^0} + \overline{u' \nabla u'} = \nabla \bar{p} + \bar{f} ; \nabla \cdot u^0 = 0 \text{ in } \Omega \times ]0,T[ \\
u^0(x,0) = u_o ; u|_{\Gamma} = u_l .
\end{cases}
$$

In the eddy viscosity models one assumes that the Reynolds stresses $\overline{u^0 \nabla u'}$ + $\overline{u' \nabla u^0}$ + $\overline{u' \nabla u'}$ are equivalent to some non linear operator of the form
$\nabla \cdot [A(\nabla u, \frac{\partial}{\partial t} \nabla u)\nabla u]$ where A is an appropriate tensor, and the main difficulty
with this approach is to choose A properly.

However we observe that while the Reynolds stresses are difficult to model,
the homogeneous part of the turbulence, u', is not difficult to model. Indeed,
from KOLMOGOROV (1941) first similarity hypothesis u' can be modelled from a
universal random process, w, scaled by the viscosity v and local rate of
dissipated energy v|∇u°|² of the flow we wish to simulate, filtered to
remove the eddies and frequencies which are not in the inertial subrange
of the turbulence considered, and transported by the mean flow at speed u°.
Therefore, as shown in § 1, if ε⁻¹ denotes the largest wave number of the
inertial subrange and ε² the largest allowed frequency, u' is found to be
in the form

$$
(4) \quad u'(x,t) = v^\epsilon (v, |\nabla u^0|^2 (x,t), \frac{x-u^0 t}{\epsilon}, \frac{t}{\epsilon^2})
$$

where $v^\epsilon$ is a highly oscillating random process.

Then our problem becomes : find the mean value <u<sup>&epsilon;</sup>> of u<sup>&epsilon;</sup> solution of

$$
(5) \quad \left\{
\begin{array}{l}
\displaystyle \frac{\partial u}{\partial t} - v \Delta u^{\epsilon} + (u^{\epsilon} + v^{\epsilon}) \nabla (u^{\epsilon} + v^{\epsilon}) = \nabla p^{\epsilon} + \bar{f} - \frac{\partial v^{\epsilon}}{\partial t} + v \Delta v^{\epsilon} ; \nabla \cdot u^{\epsilon} = 0 \\[1em]
\displaystyle u^{\epsilon}(t=0) = u^{0}, \quad u^{\epsilon}|_{\Gamma} = u_{1}
\end{array}
\right.
$$
where $v^\epsilon$ is a given stochastic process. If $v^\epsilon$ simulates $u'$ correctly then $\langle u^\epsilon \rangle$ should be close to $u^0$. Therefore, we can expect that $u^\epsilon = u^0 + \epsilon u^1 + \epsilon^2 u^2 + \dots$ so that our main concern will be to study the limit of $u^\epsilon$ when $\epsilon$ tends to zero.

From the numerical point of view (5) is still very difficult because $v^\epsilon$ is a highly oscillating function of $x$ and $t$ when $\epsilon \ll 1$ and many points of discretization will be needed as for (1). Unlike (3), it is useless to average (5) because it would yield a still more difficult integro-differential equation. So we must turn to some more sophisticated method of averaging.

The first method, from the theory of stochastic differential equation, would be to find the law of probability of $u^\epsilon$, which strictly speaking is also a random process. In §2 we discuss this point of view and the very interesting work of STRATONOVITCH (1963), PAPANICOLAU-KOLHER (1974) and the references therein, and of MEYERS-O'BRIEN-SCOTT (1977) as applied to the averaging of the equation for the temperature $\theta^\epsilon$:

$$ (6) \quad \begin{cases} \frac{\partial \theta^\epsilon}{\partial t} + (u^0 + v^\epsilon) \nabla \theta^\epsilon = 0 \\ \theta^\epsilon(t=0) = \theta_o, \quad \theta_n^\epsilon|_T = \theta_l \end{cases} $$

which, when $v^\epsilon$ is an oscillating random process, offers the same difficulties as (5). But in the above references it is shown (in a different context) that when

$$ (7) \quad v^\epsilon = \frac{1}{\epsilon} v(x,t, \frac{t}{\epsilon^2}) ; \langle v \rangle = 0, $$

then $\langle \theta^\epsilon \rangle$ tends to $\langle \theta^0 \rangle$ solution of

$$ (8) \quad \frac{\partial \langle \theta^0 \rangle}{\partial t} + u^0 \nabla \langle \theta^0 \rangle + \nabla \cdot (A \nabla \langle \theta^0 \rangle) = 0 ; \langle \theta^0 \rangle (t=0) = \theta_o ; \langle \theta^0 \rangle |_T = \theta_l, $$

where A is a tensor with coefficients $a_{ij}(x,t)$ computed by averaging a time integral of $v_i v_j$. This result is interesting because it supports the concept of eddy thermal conductivity. However the result rely on assumption (7) which, by the way, does not mean that $v^\epsilon$ is frequently of large amplitude, but rather that the turbulence contains very high frequencies. In §2 we discuss also the case
$$ (9) \qquad v^{\epsilon}(x,t) = \frac{1}{\epsilon} v(x,t, \frac{x-u^0 t}{\epsilon}, \frac{t}{\epsilon^2}) $$

and show formally that a result similar to (8) is to be excepted. Indeed when $v$ is periodic with respect to $\frac{x}{\epsilon}$ and $\frac{t}{\epsilon^2}$, (and $v$ is not convected by $u^0$), a result like (8) can be shown by the theory of homogenization. This brings us to the second method to average (6).

The theory of homogenization was developed in structural mechanics to study composite materials such as alloys, plastics, ... which are inhomogeneous but have a periodic structure. It shows that if the structure is very fine, the composite material is then equivalent to an homogeneous material whose characteristics can be computed. The theory is more general and applies to study the limit of solutions of partial differential equation with highly oscillating coefficient (see BENSOUSSAN-LIONS-PAPANICOLAU (1978) and the bibliography therein).

In §3 we apply the theory of homogenization to (5) and show that if $v^\epsilon$ is a deterministic function given by (9) and periodic with respect to $(x-u^0t)/\epsilon$ and $t/\epsilon^2$ on $C = \mathbb{Z}x]0,S[ \text{ then } u^\epsilon$ solution of (5) tends to

$$ (10) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial u^0}{\partial t} - v \Delta u^0 + u^0 \nabla u^0 + B(u^0, \nabla u^0, \nabla(\nabla u^0)) = f + \nabla p, \quad \nabla \cdot u^0 = 0 \\[1.5ex] \displaystyle u^0(t=0) = u^0; \quad u^0|_{T} = u^1, \end{array} \right. $$

where $B$ is a linear function of $u^0, \nabla u^0, \nabla(\nabla u^0)$ and can be computed as the mean on the cell of periodicity $C$ of $v$ of $v\nabla u^1 + u^1\nabla v$ where $u^1$ is solution on $C$ of

$$ (11) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial u^1}{\partial \tau} - v \Delta_z u^1 + u^1 \nabla_z v + v \nabla u^1 = \nabla q - v \nabla_x u^0 - u^0 \nabla_x v; \quad \nabla_z \cdot u^1 = 0, \quad \text{on } C \\[1.5ex] \displaystyle u^1 \quad \text{periodic on } C. \end{array} \right. $$

If $v$ is a function of the local rate of dissipated energy then $B$ is a non linear function of $\nabla u^0$.
In (11) z stands for $(x-u^0t)/\epsilon$ and $\tau$ for $t/\epsilon^2$. This result is analogous to
(8) because when $v$ tends to zero and (11) is integrated on its characteristics then $v\nabla u^1 + u^1\nabla v$ involves integrals of $v_i v_j$. Note also that these results seem to indicate that if all wave numbers smaller than $\epsilon^{-1}$ are filter to construct $v^\epsilon$ then all frequencies smaller than $\epsilon^{-2}$ should be filtered also, otherwise an equation like (11) cannot be obtained.

The result above is derived formally only because it is very difficult to handle
mathematically the transport variable $x-u^0t$. Although (10) is the full result,
because it is difficult to handle, we felt like developping a simpler model
where the transport/interaction term $u^\epsilon \nabla v^\epsilon$ in (5) is neglected. Therefore
in §4 we restrict ourselves to the case where

$$
(12) \quad v^{\epsilon}(x,t) = \frac{1}{\epsilon} v(x,t, \frac{x}{\epsilon}, \frac{t}{\epsilon^2}), \quad v \text{ periodic with respect to } (y,\tau) = \left(\frac{x}{\epsilon}, \frac{t}{\epsilon^2}\right) \text{ on C.}
$$

and study the limit $u^0$ of $u^\epsilon$ solution of

$$
(13) \quad \begin{cases}
\displaystyle \frac{\partial u^\epsilon}{\partial t} - v \Delta u^\epsilon + (u^\epsilon + v^\epsilon) \nabla u^\epsilon = f + \nabla p, & \nabla \cdot u^\epsilon = 0 \text{ in } \Omega \times ]0, T[ \\
u^\epsilon(t=0) = 0, \quad u^\epsilon|_{T} = 0.
\end{cases}
$$

Then we prove that if $v^\epsilon$ is periodic with respect to the last variables, $u^\epsilon$
tends to $u^0$ solution of

$$
(14) \quad \left\{
\begin{array}{l}
\displaystyle \frac{\partial u^0}{\partial t} - v \Delta u^0 + u^0 \nabla u^0 + \nabla \cdot (A \nabla u^0) = \nabla p + f \quad \nabla \cdot u^0 = 0 \quad \text{in } \Omega \times ]0, T[ \\
u^0(t=0) = 0, \quad u^0|_{T} = 0
\end{array}
\right.
$$

where $-A$ is a positive definite tensor, function of $(x,t)$, and therefore possibly
of $v|\nabla u^0|^2$, such that $A\nabla u^0$ is the mean of $v \otimes u^1$ on $C$ and $u^1$ is the solution of

$$
(15) \quad \left\{
\begin{array}{l}
\displaystyle \frac{\partial u^1}{\partial \tau} - v \Delta_y u^1 + v \nabla_y u^1 = \nabla q + v \nabla_x u^0, \quad \nabla_y \cdot u^1 = 0 \quad \text{in } C \\
u^1 \quad \text{periodic in } C.
\end{array}
\right.
$$

Then we conclude in §5 with some comments on the possible experimental and
numerical check of these results. These results seem reasonable but they
are based on the peculiar behavior of $v^\epsilon$, in (9), with respect to $\epsilon$. In §1
we show that the part of the inertial subrange which is just behind the
cut-off wave number $\epsilon^{-1}$ is indeed of the form

$$ (16) \qquad v^{\epsilon}(x,t) = \epsilon^{\beta} v(x,t, \frac{x-u^0 t}{\epsilon}, \frac{t}{\epsilon}) $$

for some $\beta$ (depending on the hypothesis on the behavior of $v^\epsilon$ with respect to $t$), where $v$ is a random process periodic with respect to the last variables. Therefore the best way to check the present point of view would be to measure experimentally (or simulate it) an homogeneous turbulence, scale it and filter it to obtain a set of values to generate the random process $v^\epsilon$ then check the results on the temperature equation by computing numerically $\theta^0$ solution of

$$ \frac{\partial \theta^0}{\partial t} + u^0 \nabla \theta^0 + \lambda^\epsilon \nabla \cdot (A \nabla \theta^0) = e ; \quad \theta^0(t=0) = \theta_0, \theta^0|_T = \theta_1 $$

for different values of $\lambda^\epsilon$ (to be adjusted) and where $A$ is a random matrix whose coefficients are the means of $v_i v_j$ integrated on the v-characteristics. Then a similar check can be made on $u^0$ by (10) or (14). For $u^0$ it is more convenient to keep $A$ random in order to destabilize numerically the unstable large structures.
# 1. - MODELLING OF THE HOMOGENEOUS PART OF THE TURBULENCE

The Kolmogorov similarity hypothesis assumes that if a turbulence is locally isotropic and homogeneous it can be expressed in terms of the viscosity $\nu$, the local rate of energy dissipated $e$ and a universal random field $w$. More precisely (see MONIN and YAGLOM (1975, p. 348)) the velocity of the flow is given by

$$ (1.1) \quad u(x_o+x)_{t_o+t} = u(x_o,t_o) + (\nu e)^{1/4} w(x(\frac{e}{\nu}), t(\frac{e}{\nu})^{1/2}) $$

where

$$ (1.2) \quad e = \frac{\nu}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)^{1/2} $$

and $\bar{e}$ denotes the mean value.

In general only a part of $u$ is an isotropic homogeneous turbulence, namely the eddies of the "inertial subrange" of the spectrum of $u$. Furthermore this turbulence is convected by the mean flow.

Therefore, in a frame of reference moving at the mean velocity of the flow $u^o$, we must filter all the eddies whose wave number is smaller than the lower bound, say $\epsilon^{-1}$, of the inertial subrange.

This operation can be performed by a convolution of $w$ with the inverse Fourier transform $G(\frac{x}{\epsilon})$ of the step function $Y(\epsilon\vec{k})$:

$$ (1.3) \quad Y(\epsilon\vec{k}) = \begin{cases} 1 & \text{if } |\vec{k}| < \epsilon^{-1} \\ G(\frac{x}{\epsilon}) = \prod_{i=1}^{3} \frac{\sin(\frac{x_i}{\epsilon})}{\frac{x_i}{\epsilon}} & \\ 0 & \text{otherwise.} \end{cases} $$

We shall see later that if a filter in $\frac{x}{\epsilon}$ is used it is necessary to filter in time also the frequencies larger than $\epsilon^{-2}$. Therefore we shall use

$$ (1.4) \quad G(\frac{x}{\epsilon}, \frac{t}{\epsilon^2}) = \frac{\sin(\frac{t}{\epsilon^2})}{\frac{t}{\epsilon^2}} G(\frac{x}{\epsilon}) . $$

Then the homogeneous part of the turbulence that we will use is
$$ (1.5) \qquad v^\epsilon(x,t) = (\nu e)^{1/4} \left[ w\left(x' \left(\frac{e}{\nu}\right)^{1/4}\right), t\left(\frac{e}{\nu}\right)^{1/2} \right] - w\left(x' \left(\frac{e}{\nu}\right)^{1/4}\right), t\left(\frac{e}{\nu}\right)^{1/2} \left( G\left(\frac{x'}{\epsilon}, \frac{t}{\epsilon^2}\right) \right) $$

where is the operator of convolution and

$$ (1.6) \qquad x' = x - u^0 t. $$

Note that we have neglected the dissipative subrange or rather we have assumed that it is contained in w.

For practical purposes $v^\epsilon$ can be generated by recording w on tapes from an experiment (say a flow behind a grid), scaling and filtering.

From (1.5) the dependence of $v^\epsilon$ with respect to $\epsilon$ is not clear. In the following paragraphs we shall assume that $v^\epsilon$ depends of $\epsilon$ as

$$ (1.7) \qquad v^\epsilon(x,t) = \frac{1}{\epsilon} v(x,t, \frac{x-u^0 t}{\epsilon}, \frac{t}{\epsilon^2}) $$

where v is a zero mean stationary random process or a deterministic function periodic with respect to the last variables $(\frac{x^2}{\epsilon}, \frac{t}{\epsilon^2})$ on some cell $C = Y \times [0, S]$, and with zero mean value on C. The coefficient $\frac{1}{\epsilon}$ in (1.7) is strange indeed and we can produce only one evidence to justify its use, beside the fact that the results we shall derive are physically meaningful.

The argument rely on the interesting following property : if $\eta(t)$ is any stationary zero mean random process then, when $\epsilon$ tends to zero, $\frac{1}{\epsilon} \eta(\frac{t}{\epsilon^2})$ converges to a zero mean Gaussian white noise, by the central limit theorem. Therefore (1.7) does not mean that very large fluctuations are frequent, it means that high frequencies are frequent, because a white noise is characterized by the lack of correlations between its values at two different times.

In §3 and 4 we shall also assume that v is periodic with respect to $(x-u^0 t)/\epsilon$ and eventually with respect to $t/\epsilon^2$.

Such hypothesis is not in contradiction with Kolmogorov's second 1941 - similarity hypothesis which implies that
$$
(1.8) \quad v^{\varepsilon}(x,t) = \varepsilon^{\beta} v(x,t, \frac{x-u^0 t}{\varepsilon}) * \left[ \frac{\sin(t/\varepsilon^2)}{t/\varepsilon^2} \right]
$$

for some $\beta$ and $v$, and that the part of the spectrum of $v^\epsilon$ just behind the cut-off
wave number $1/\epsilon$ is periodic with respect to the last variables. Note that if
we knew the behaviour of turbulence with respect to time (Taylor's hypothesis ?)
we could, by a similar argument, show that $v^\epsilon$ depend of $\epsilon$ through $\frac{x-u_0 t}{\epsilon}$, $\frac{t}{\epsilon^2}$
and a multiplicative factor $\epsilon^\beta$.

To show (1.8) we denote by $\hat{w}$ the Fourier transform of $w$ with respect to $x$,
and since we work at a fixed time $t$ we omit to filter in time and omit to
write the parameter $t$. Then the inertial subrange is also (provided the stochastic
integral exists)

$$
(1.9) \quad v^{\epsilon}(x) = \int_{\left|k\right| \ge \frac{1}{\epsilon}} \exp\left(z, \pi k x'\right) d\hat{w}(k, x)
$$

where $x' = x - u_0 t$. From the second similarity hypothesis of Kolmogorov, the second
order moments of $\hat{w}$ are power functions of $k$, therefore in the stochastic integral
(1.9) $d\hat{w}(\epsilon^{-1}k, x)$ may be replaced by $\epsilon^\beta d\hat{w}(k, x)$ for some $\beta$. Hence letting $k' = \epsilon k$

$$
(1.10) \quad v^{\epsilon}(x) = \epsilon^{\beta} \int_{|k'| \ge 1} \exp\left(2i\pi k' \frac{x}{\epsilon}\right) d\hat{w}(k', x) .
$$

Now restrict the domain of integration in (1.10) to $1+\epsilon M \ge |k'| \ge 1$, then the
integral in (1.10) is taken over a narrow strip and can be approximated by $\epsilon$
time the boundary integral :

$$
(1.11) \quad v_c^\epsilon(x) = \epsilon^{\beta+1} \int_{1+\epsilon M \ge |k'| \ge 1} \exp\left(2i\pi k' \frac{x}{\epsilon}\right) d\hat{w}(k',x) = \epsilon^{\beta+1} \int_{|k'|=1} \exp\left(2i\pi k' \frac{x}{\epsilon}\right) d\hat{w}(k',x)
$$

Finally if the norm chosen is $|k'| = \sup_{i=1,2,3} |k'_i|$, $v_c^\epsilon$ is a linear combination
of periodic functions :

$$
v_c^\epsilon(x) = \epsilon^{\beta+1} \sum_j \exp\left(2i\pi \left(\frac{x_j}{\epsilon}\right)^{\frac{1}{\beta}}\right) \int_{k'' \in S_j} \exp\left(2i\pi k'' \frac{x}{\epsilon}\right) d\hat{w}\left(\ell^{j} + k'', x\right)
$$

where $S_j$ are all the faces of the cube $|k'| = 1$ and $\ell^j$ the corresponding normal
unit vector to $S_j$. Letting again $k_l = k''\epsilon^{-1}$ and assuming that $\hat{w}$ is stationary
in $k_l$ we find
$$ (1.12) \quad v_c^\epsilon(x) \approx \epsilon \sum_j \exp(2i\pi \vec{l}^j, \frac{\vec{x}'}{\epsilon}) \int_{k_1 \in (\vec{l}^j)^\perp} \exp(2i\pi \vec{k}, \vec{x}') d\hat{w}(k, x) $$

Of course this derivation is purely formal.
## 2. AVERAGING OF THE TEMPERATURE EQUATION

For incompressible flows in standard conditions (LANDAU-LIFCHITZ (1953, Eq. 50,2)), the temperature $θ(x,t)$ is the solution of

$$ (2.1) \qquad \begin{cases} \frac{\partial \theta}{\partial t} - k \Delta \theta + u \nabla \theta = e \quad \forall (x,t) \in \Omega \times ]0,T[ = Q \\ \theta(x,0) = \theta_0(x) \quad \forall x \in \Omega \quad ; \quad \theta(x,t) = \theta_1(x,t) \quad \forall (x,t) \in \Gamma \times ]0,T[ \end{cases} $$

where $k$ is the thermal conductivity ($k<<1$), $u(x,t)$ is the velocity of the flow, $e(x,t)$ the rate of viscous energy in the flow; $\Omega$ is the region occupied by the fluid, $\Gamma$ is its boundary and $(\theta_0, \theta_1)$ is a set of boundary conditions.

For the time being we shall ignore that the homogeneous part of the turbulence is convected by the mean flow. Then according to §1 we choose

$$ (2.2) \qquad u(x,t) = u^o(x,t) + \frac{1}{\epsilon} v(x,t, \frac{x}{\epsilon}, \frac{t}{\epsilon^2}) + O(\epsilon) $$

where $u^o$ is the mean flow and $\frac{1}{\epsilon}v$ is the inertial subrange; $O(\epsilon)$ represents the coupling between $u^o$ and $\frac{1}{\epsilon}v$ which as we shall show in §3 is indeed such that

$$ (2.3) \qquad |O(\epsilon)| \le K\epsilon. $$

The problem we wish to solve is:

1.°) if $v$ is a zero mean stationary random process find

$$ (2.4) \qquad \langle \theta^o \rangle (x,t) = \lim_{\epsilon \to 0} \langle \theta^\epsilon \rangle $$

where $\langle \theta^\epsilon \rangle$ is the expected value of the solution of (2.1) with $u$ given by (2.2).

2.°) if $v$ is a periodic deterministic function with zero mean value on the cell of periodicity, find

$$ (2.5) \qquad \theta^o(x,t) = \lim_{\epsilon \to 0} \theta^\epsilon. $$
To simplify the computations we shall assume that

$$ (2.6) \qquad \nabla \cdot v^{\epsilon} = 0 \quad \forall \epsilon \text{ i.e. } \nabla_X \cdot v(x,t,y,\tau) = 0, \nabla_Y \cdot v(x,t,y,\tau) = 0 \quad \forall x,t,y,\tau. $$

To begin, we recall an interesting result of MEYER-O'BRIEN-SCOTT (1977) : if $k=0$ and if the random process $v$ does not depend upon $\frac{x}{\epsilon}$ then, provided that the "mixing" hypothesis of PAPANICOLAU-KOHLER (1974, (2.1)-(2.7)) holds, the mean $\langle \theta^{\epsilon} \rangle$ converges to $\langle \theta^0 \rangle$ solution of

$$ (2.7) \qquad \begin{cases} \frac{\partial \langle \theta^0 \rangle}{\partial t} + u^0 \nabla \langle \theta^0 \rangle + \nabla \cdot (A \nabla \langle \theta^0 \rangle) = e \text{ in } Q \\ \langle \theta^0 \rangle(x,0) = \theta_0(x) \quad \forall x \in \Omega \quad ; \quad \langle \theta^0 \rangle|_{\Gamma} = \theta_1, \end{cases} $$

where

$$ (2.8) \qquad a_{ij}(x,t) = \lim_{\epsilon \to 0} \epsilon^{-3} \int_{t}^{t+\epsilon} \int_{t}^{t_2} \langle v_i(x,t, \frac{t_1}{\epsilon}) v_j(x,t, \frac{t_2}{\epsilon}) \rangle dt_1 dt_2 $$

where $\langle \rangle$ is the expected value with respect to the probability law of $v$.

The proof of (2.7)-(2.8) is based upon an integration of (2.1) by the method of characteristics. One study the limit of the solution of

$$ (2.9) \qquad \frac{dx}{dt} = u^0(x^\epsilon, t) + \frac{1}{\epsilon} v(x^\epsilon, t, \frac{t}{\epsilon}) $$

As shown in PAPANICOLAU-KOHLER (1974) (see also STRATONOVITCH (1963)) $x^\epsilon$ converges to a diffusion process and therefore $\langle \theta^0 \rangle$ can be found from Itô's calculus to satisfy (2.7).

This result justifies the concept of turbulent thermal conductivity.

Note that (2.7) is easy to simulate numerically (by opposition to (2.1)) because $A$ is a positive definite well conditioned matrix and there are no oscillatory coefficients.

Now we wish to extend this result to the case where $v$ depends upon $x/\epsilon$ also.
We begin with a formal argument : assume that v is deterministic and that $\theta^\epsilon$
can be expanded as

$$
(2.10) \quad \theta^{\epsilon}(x,t) = \theta^{0}(x,t) + \epsilon\theta^{1}(x,t,\frac{x}{\epsilon},\frac{t}{\epsilon}) + \epsilon^{2}\theta^{2}(x,t,\frac{x}{\epsilon},\frac{t}{\epsilon}) + \dots
$$

then (2.1) implies

$$
(2.11) \quad \left[ \left( \frac{\partial}{\partial t} + \epsilon^{-2} \frac{\partial}{\partial \tau} \right) - k \left( \epsilon^{-2} \Delta_y + \epsilon^{-1} \Delta_{xy} + \Delta_x \right) + (u^0 + \epsilon^{-1} v) \left( \epsilon^{-1} \nabla_y + \nabla_x \right) \right] \cdot \left( \theta^0 + \epsilon \theta^1 + \epsilon^2 \theta^2 + \dots \right) = 0
$$

where

$$
(2.12) \quad y = \frac{x}{\epsilon} ; \tau = \frac{t}{\epsilon^2} ; \Delta_{xy} = 2 \frac{\partial^2}{\partial x_i \partial y_i} .
$$

If we write that (2.11) holds for all (x,t,y,τ,ε) instead of all (x,t,ε),
we find that

$$
(2.13) \quad \frac{\partial \theta^1}{\partial \tau} - k \Delta_y \theta^1 + v \nabla_y \theta^1 + v \nabla_x \theta^0 = 0
$$

$$
(2.14) \quad \frac{\partial \theta^0}{\partial t} - k \Delta_x \theta^0 + u^0 \nabla_x \theta^0 + v \nabla_x \theta^1 + u^0 \nabla_y \theta^1 - k \Delta_{xy} \theta^1 + \frac{\partial \theta^2}{\partial \tau} - k \Delta_y \theta^2 + v \nabla_y \theta^2 = e
$$

If $k=0$ (2.13) can be integrated on the $(y,\tau)$ characteristics

$$
(2.15) \quad \frac{dy}{d\tau} = v(x,t,y,\tau)
$$

$$
(2.16) \quad \theta^1(x,t,y_o + \int_o^\tau y(\sigma)d\sigma, \tau) = - \left[ \int_o^\tau v(x,t,y(\sigma),\sigma)d\sigma \right] \nabla_x \theta^0 .
$$

Loosely speaking, v is an oscillatory function with zero mean value, therefore
$\theta^1$ is also oscillatory with zero mean value, similarly (2.14) implies that $\theta^2$
has zero mean value. Therefore if we average (2.14) with respect to $(y,\tau)$ and
denote by $-y\tau$ the averaging operation, we find that

$$
(2.17) \quad \frac{\partial \theta^0}{\partial t} + u^0 \nabla_x \theta^0 + v \nabla_x \theta^1 = e
$$

and from (2.16)

$$
(2.18) \quad \overline{v \nabla_x \theta^1} = \frac{\partial}{\partial x_i} \overline{v_i \theta^1} = - \frac{\partial}{\partial x_i} \left[ \left( v_i \int_{(\tau, y(\tau))}^{y\tau} v_j d\sigma \right) \frac{\partial \theta^0}{\partial x_j} \right] .
$$
So we find that the results (2.7)(2.8) seems also valid in the general case
(but k=0) provided that the time integration in (2.8) be replaced by an
integration on the stream lines of v and the average be taken over (y,τ).
... seems also valid in the general case

Of course all this would be difficult to prove rigourously. However if v is
random and periodic with respect to y on a cell Y = ]0,Y₁[×]0,Y₂[×]0,Y₃[ it
seems that (2.8) or a similar result can be shown ; the proof is in progress.
If v is deterministic and periodic with respect to (y,τ) on some cell
C = Y × ]0,S[, Y = ]0,Y₁[×]0,Y₂[×]0,Y₃[ i.e. if

$$
(2.19) \quad v(x,t,y_i+Y_i, \tau+S) = v(x,t,y_i, \tau) \quad \forall (x,t,y,\tau) \in \Omega \times ]0,T[ \times C
$$

then the following results hold

Theorem 1 : If k > 0, v satisfies (2.19) and ∇_x·v = ∇_y·v = 0; If

$$
\int_Y v \, dy = 0 \quad \forall (x,t,\tau) \in Q \times ]0,S[ ,
$$

and v=0 in a neighborhood of Γ and t=0, and if v ∈ L^∞(Q×C), then the solution
θ^ε of (2.1) with u given by

$$
(2.20) \quad u^{\epsilon}(x,t) = u^{0}(x,t) + \frac{1}{\epsilon} v \left( x, t, \frac{x}{\epsilon}, \frac{t}{\epsilon^2} \right)
$$

converges to θ° in the sense

$$
(2.21) \quad \left( \int_{\Omega} |\theta^{\epsilon} - \theta^{0}|^2 dx \right)^{1/2} \le \frac{K}{\kappa} \epsilon, \left( \int_{\Omega} |\nabla(\theta^{\epsilon} - \theta^{0} - \epsilon\theta^{1})|^2 dx \right)^{1/2} \le \frac{K}{\kappa} \epsilon, \forall \epsilon, \forall t \in ]0, T[
$$

where θ°,θ¹ are the solutions of

$$
(2.22) \quad \begin{cases} \frac{\partial \theta^0}{\partial t} - k \Delta \theta^0 + u^0 \nabla \theta^0 + \nabla \cdot (A \nabla \theta^0) = e \text{ in } \Omega \times ]0, T[ \\ \\ \theta^0(x,o) = \theta_o(x) \text{ in } \Omega ; \quad \theta^0|_{\Gamma} = \theta_I . \end{cases}
$$

$$
(2.23) \quad \theta^1(x,t,y,\tau) = \bar{\Sigma}(x,t,y,\tau) \nabla \theta^0(x,t)
$$

where Σ is the solution of
Now define

$$
(2.27) \quad \phi^{\epsilon}(x,t) = \theta^{\epsilon}(x,t) - \theta^{0}(x,t) - \epsilon \theta^{1}(x,t, \frac{x}{\epsilon}, \frac{t}{\epsilon}) - \epsilon^{2} \theta^{2}(x,t, \frac{x}{\epsilon}, \frac{t}{\epsilon})
$$

where $\theta^2$ is the solution of

$$
(2.28) \quad \left\{
\begin{array}{@{}l@{}}
  \displaystyle \frac{\partial \theta^2}{\partial t} - k \Delta_y \theta^2 + v \nabla_y \theta^2 = - u^o \nabla_y \theta^1 - v \nabla_x \theta^1 + \overline{v \nabla_x \theta^1}^{1y\tau} + k \Delta_{xy} \theta^1 \\[2ex]
  \displaystyle \theta^2(y, \tau) \quad \text{periodic on } C \; ; \quad \overline{\theta^2}^{y\tau} = 0.
\end{array}
\right.
$$

Then putting (2.77) into (2.1) we find that

$$
(2.29) \quad \left\{
\begin{array}{@{}l@{}}
  \displaystyle \left. \frac{\partial \phi^\epsilon}{\partial t} - k \Delta \phi^\epsilon + \left(u^o + \frac{1}{\epsilon} v\right) \nabla \phi^\epsilon \right|_{\substack{y=x/\epsilon \\ \tau=t/\epsilon^2}} = (\epsilon k \Delta_x \theta^1 + \epsilon k \Delta_{xy} \theta^2 + \epsilon^2 k \Delta_x \theta^2 - \right. \\
  \displaystyle \quad \left. - \epsilon^2 \frac{\partial \theta^2}{\partial t} - \epsilon u^o \nabla_x \theta^1 - \epsilon^2 u^o \nabla_x \theta^2 - \epsilon v \nabla_x \theta^2 - \epsilon u^o \nabla_y \theta^2) \right|_{\substack{y=x/\epsilon \\ \tau=t/\epsilon^2}} \\[2ex]
  \phi^\epsilon(x,o) = 0, \quad \phi^\epsilon|_T = 0
\end{array}
\right.
$$

Multiplying (2.29) by $\phi^\epsilon$ and integrating over $Q$ yields

$$
(2.30) \quad \frac{1}{2} \int_{\Omega} \phi^{\epsilon} (x,t)^2 dx + k \int_{\Omega} |\nabla \phi^{\epsilon}|^2 (x,t) dx \le \epsilon K \left( \int_Q |\nabla \phi^{\epsilon}|^2 dx dt \right)^{1/2}
$$

where $K$ is the $L^2(0,T;H^{-1}(\Omega))$-norm of

$$
\left[ k \Delta_x \theta^1 + \epsilon (k \Delta_{xy} \theta^2 - u^o \nabla_x \theta^1 - v \nabla_x \theta^2 - u^o \nabla_y \theta^2) + \epsilon^2 (k \Delta_x \theta^2 - \frac{\partial \theta^2}{\partial t} - u^o \nabla_x \theta^2) \right] \Big|_{\substack{y=\frac{x}{\epsilon t} \\ \tau=\frac{\epsilon^2}{\epsilon t}}} = 0
$$

It is now straightforward to deduce (2.21) from (2.30).

**Remark 2.4** : To show that theorem 1 holds also when k=0 one must show that (2.24) has a periodic solution (this puts some restrictions on v) and that the solution is regular.
**Remark 2.5** : If v is zero on the boundary of Y and k=0 the coefficients $a_{ij}$ of (2.25) can be computed exactly. Indeed by multiplying (2.24) by $y_i$ and integrating over C we find that

$$ (2.31) \quad \overline{v_i \cdot \Sigma_j}^{-y\tau} = \overline{v_j \cdot y_i}^{-y\tau}. $$

Therefore when $\nabla_x \cdot v = 0$

$$ (2.32) \quad \nabla \cdot A \nabla \theta^0 = \frac{\partial}{\partial x_i} \left( \overline{v_j \cdot y_i} \frac{\partial \theta^0}{\partial x_j} \right) = \sum_{j \ge i} \frac{\partial^2}{\partial x_i \partial x_j} \left[ \overline{(v_j \cdot y_i) \theta^0} \right]. $$

But

$$ (2.33) \quad \overline{v_j \cdot y_i}^{-y\tau} + \overline{v_i \cdot y_j}^{-y\tau} = v \nabla (y_i y_j) = 0 $$

therefore in that case

$$ (2.34) \quad \nabla \cdot (A \nabla \theta^0) = 0. $$

This result is not surprising because when $v|_{\partial Y} = 0$, the stream lines are closed and therefore from the divergence theorem (2.18) gives also $v \nabla_x \theta^1 = 0$.

However it is interesting to compare (2.32) with the crude average of (2.1) based on the hypothesis that $\nabla \theta$ varies according to

$$ \nabla \theta(x+\epsilon y) \approx \nabla \theta(x) + \epsilon y \nabla(\nabla \theta) $$

then if $B(x, \epsilon)$ denotes the sphere of radius $\epsilon$ and center x

$$ (2.35) \quad \int_{B(x, \epsilon)} u(z) \nabla \theta(z) dz \approx u \nabla \theta + \left(\frac{\epsilon}{4\pi} \int_{B(0,1)} u_i(x+\epsilon y) y_j dy\right) \frac{\partial^2 \theta}{\partial x_i \partial x_j}. $$
### 3. AVERAGING OF THE NAVIER STOKES EQUATION

To simulate numerically the solution of

$$ (3.1) \quad \begin{cases} \displaystyle \frac{\partial u}{\partial t} - v \Delta u + u \nabla u = \nabla p + f ; \nabla \cdot u = 0 \text{ in } \Omega \times ]0,T[ = Q \\ u(t=0) = u_o \text{ in } \Omega ; u|_{\Gamma} = u_l \quad t \in ]0,T[ \end{cases} $$

we assume that the untractable part of $u$ by the computer is a stochastic process $v^\epsilon$; more precisely we assume that the solution of (3.1) is of the form

$$ (3.2) \quad u(x,t) = u^\epsilon(x,t) + \frac{1}{\epsilon} v(x,t, \frac{x-u^0 t}{\epsilon}, \frac{t}{\epsilon^2}) $$

where $u^0$ is the mean flow and is related to $u^\epsilon$ by

$$ (3.3) \quad \lim_{\epsilon \to 0} u^\epsilon(x,t) = u^0(x,t) \quad (x,t) \in Q. $$

The random process $v$ may also depend of the local rate of dissipated energy $e = v|\nabla u^0|$.

As we pointed out in the introduction, the study of the limit of $u^\epsilon$ when $\epsilon \to 0$ and when $v$ is a random process is much too difficult. Therefore we shall assume instead that $v$ is a deterministic periodic function in order to derive some PDE for $u^0$ and as we have shown in §2 it is not absurd to conjecture that the result holds in the general case.

#### Assumptions

(3.4) $v(x,t,z,\tau)$ is periodic in $(z,\tau)$ on $Z \times ]0,S[$ = C

where $Z = ]0,Z_1[\times ]0,Z_2[\times ]0,Z_3[$

$$ (3.5) \quad \int_Z v \, dz = 0 \quad \forall (x,t,\tau) \in Q \times ]0,S[$$

$$ (3.6) \quad \nabla_z \cdot v = 0 $$

(3.7) $v=0$ in the neighborhood of the boundaries $\Gamma$ and at $t=0$.
Notice that we do not assume that $\nabla_x \cdot v = 0$ because it may be difficult to model $v$ by (1.5) while satisfying $\nabla_x \cdot v = 0$.

Letting

$$ (3.8) \quad v^{\epsilon}(x,t) = \frac{1}{\epsilon} v(x,t, \frac{x-u_0^t}{\epsilon}, t/\epsilon^2) $$

and replacing $u$ given by (3.2) in (3.1), we find

$$ (3.9) \quad \left\{ \begin{array}{l} \displaystyle \frac{\partial u^{\epsilon}}{\partial t} - v \Delta u^{\epsilon} + (u^{\epsilon} + v^{\epsilon}) \nabla u^{\epsilon} + u^{\epsilon} \nabla v^{\epsilon} = \nabla p + f - \left(\frac{\partial v^{\epsilon}}{\partial t} - v \Delta v^{\epsilon} + v^{\epsilon} \nabla v^{\epsilon}\right) \\[1.5ex] \displaystyle \nabla \cdot u^{\epsilon} = - \nabla \cdot v^{\epsilon} \quad ; \quad u^{\epsilon}(t=0) = u_o, \quad u^{\epsilon}|_{\Gamma} = u_1. \end{array} \right. $$

We may assume that (3.9) is the equation that defines $u^\epsilon$ once $v^\epsilon$ is known. Therefore our problem fits into the framework of homogenization theory which studies the limits of solutions of PDEs, with highly oscillating periodic coefficients.

For practical purposes it is convenient to write also

$$ (3.10) \quad p = p^{\epsilon} + q^{\epsilon} $$

for a given $q^\epsilon$.

Then we search for $(u^\epsilon, p^\epsilon)$ in the form

$$ (3.11) \quad u^{\epsilon}(x,t) = (\tilde{u}^0(x,t,z,\tau) + \epsilon u^1(x,t,z,\tau) + \epsilon^2 u^2(x,t,z,\tau) + \dots) \begin{cases} z = \frac{x-u^0 t}{\epsilon} \\ \tau = \frac{t}{2} \end{cases} $$

$$ (3.12) \quad p^{\epsilon}(x,t) = \epsilon^{-1} p^{-1}(x,t,z,\tau) + p^{0}(x,t,z,\tau) + \epsilon p^{1}(x,t,z,\tau) + \dots \begin{cases} z = \frac{x-u^{0} t}{\epsilon} \\ = \frac{t}{2} \end{cases} $$

where $\tilde{u}^0, u^1, u^2, \dots, p^{-1}, p^0, p^1, \dots$ are periodic with respect to $(z, \tau)$ on $C$.
Now we replace $(u^\epsilon, p^\epsilon)$ in (3.9). However to compute the derivatives of $(u^\epsilon, p^\epsilon)$ we shall treat $u^0$ in the definition of $z$ as a constant. Indeed when $u^0$ is not constant the proper charge of variable is not $z$ but some integral of $u^0$ on the stream lines. Therefore we find for example that

$$ (3.13) \quad \nabla u^\epsilon = (\nabla_x + \frac{1}{\epsilon} \nabla_z)(\tilde{u}^0 + \epsilon u^1 + \epsilon^2 u^2 + \dots) $$

while

$$ (3.14) \quad \frac{\partial u}{\partial t} = (\epsilon^{-2} \frac{\partial}{\partial \tau} - \epsilon^{-1} u^0 \nabla_z + \frac{\partial}{\partial t})(\tilde{u}^0 + \epsilon u^1 + \epsilon^2 u^2 + \dots) $$

So (3.9) becomes

$$ \begin{cases} [\epsilon^{-2} \frac{\partial}{\partial \tau} - \epsilon^{-1} u^0 \nabla_z + \frac{\partial}{\partial t} - v(\epsilon^{-2} \Delta_z + \epsilon^{-1} \Delta_{xz} + \Delta_x)] \cdot (\tilde{u}^0 + \epsilon u^1 + \epsilon^2 u^2 + \dots) \\ + (\epsilon^{-1} v + \tilde{u}^0 + \epsilon u^1 + \epsilon^2 u^2 + \dots) (\epsilon^{-1} \nabla_z + \nabla_x) \cdot (\tilde{u}^0 + \epsilon u^1 + \epsilon^2 u^2 + \dots) \\[1.5em] (3.15) \quad \begin{aligned} &+ (\tilde{u}^0 + \epsilon u^1 + \epsilon^2 u^2 + \dots)(\epsilon^{-1} \nabla_z + \nabla_x)(\epsilon^{-1} v) = \\ &= (\epsilon^{-1} \nabla_z + \nabla_x) \cdot (\epsilon^{-1} p^{-1} + p^0 + \epsilon p^1 + q^+ + \dots) + f - [\epsilon^{-2} \frac{\partial}{\partial \tau} + \frac{\partial}{\partial t} - \epsilon^{-1} u^0 \nabla_z \\ &\quad - v(\epsilon^{-2} \Delta_y + \epsilon^{-1} \Delta_{xy} + \Delta_x)] (\epsilon^{-1} v) + \epsilon^{-2} v (\nabla_x + \epsilon^{-1} \nabla_z) v \end{aligned} \end{cases} $$

$$ (3.16) \quad (\epsilon^{-1} \nabla_z + \nabla_x) \cdot (\epsilon^{-1} v + \tilde{u}^0 + \epsilon u^1 + \epsilon^2 u^2 + \dots) = 0. $$

If we identify all the powers of $\epsilon$ and we take $q^\epsilon = \epsilon^{-2}q(s,t,z,\tau)$ we find:
terms in $\epsilon^{-3}$

$$ (3.17) \quad \frac{\partial v}{\partial \tau} - v \Delta_z v + v \nabla_z v = \nabla_y q $$

terms in $\epsilon^{-2}$

$$ (3.18) \quad \frac{\partial \tilde{u}^0}{\partial \tau} - v \Delta_y \tilde{u}^0 + v \nabla_z \tilde{u}^0 + \tilde{u}^0 \nabla_z v = \nabla_z p^{-1} + \nabla_x q + u^0 \nabla_z v - v \Delta_{xz} v + v \nabla_x v $$

terms in $\epsilon^{-1}$

$$ (3.19) \quad \begin{cases} \frac{\partial u^1}{\partial \tau} - u^0 \nabla_z \tilde{u}^0 - v \Delta_z u^1 - v \Delta_{xz} \tilde{u}^0 + v \nabla_z u^1 + \tilde{u}^0 \nabla_z \tilde{u}^0 + v \nabla_x \tilde{u}^0 \\ + \tilde{u}^0 \nabla_x v + u^1 \nabla_z v = \nabla_z p^0 + \nabla_x p^0 + \nabla_x p^{-1} - \frac{\partial v}{\partial t} - v \Delta_x v \end{cases} $$
and in (3.16)

$$
(3.20) \quad \nabla_z \cdot \tilde{u}^0 + \nabla_x \cdot v = 0
$$

. terms in ε°

$$
(3.21) \quad \left\{
\begin{aligned}
& \frac{\partial u^2}{\partial \tau} + \frac{\partial \tilde{u}^0}{\partial t} - \tilde{u}^0 \nabla_z u^1 - \nu \Delta_{xy} u^1 - \nu \Delta_x \tilde{u}^0 + \nu \nabla_z^2 u^1 \\
& \qquad + \tilde{u}^0 \nabla_x \tilde{u}^0 + u^1 \nabla_z \tilde{u}^0 + u^1 \nabla_x v + u^2 \nabla_z v = \nabla_z p^1
\end{aligned}
\right.
$$

$$
(3.22) \quad \nabla_z \cdot u^1 + \nabla_x \cdot \tilde{u}^0 = 0 .
$$

Finally the terms in ε in (2.10) give

$$
(3.23) \quad \nabla_z \cdot u^2 + \nabla_x \cdot u^1 = 0.
$$

To simplify these equations we introduce the

Notations :

$$
(3.24) \quad \bar{w}(x,t) = \frac{1}{z_1 z_2 z_3 S} \int_C w(x,t,z,\tau) dz \, d\tau
$$

$$
(3.25) \quad \bar{w}(x,t) = \frac{1}{z_1 z_2 z_3 S} \int_C w(x,t,t,\tau) dz
$$

$$
(3.26) \quad u^0(x,t) = \tilde{u}^0
$$

$$
(3.27) \quad u^{0'}(x,t) = \tilde{u}^{0} - u^{0} .
$$

With these notations (3.18) and (3.20) define u'
follows :

Definition of u�0'
(u⁰', p⁻¹) is the solution to

$$
(3.28) \quad \left\{
\begin{aligned}
& \frac{\partial u^0}{\partial \tau} - \nu \Delta_z u^0 + \nu \nabla_z u^0 + u^0 \nabla_z v = \nabla_z p^{-1} \\
& \nabla_z \cdot u^0 = - \nabla_x \cdot v \quad \text{in } C; \quad (u^0, p^{-1}) \text{ periodi}
\end{aligned}
\right.
$$
Note that (3.28) has a periodic solution only if

$$
(3.29) \quad \overline{\nabla_x q}^2 = \overline{v \nabla_x v}^2
$$

Now we note that (3.23) averaged over Z implies

$$
(3.30) \quad \nabla_x u^{\overline{l^2}} = 0 .
$$

Therefore we define $u^1$ by (3.19) and (3.22) :

Definition of $u^1$

Given $(v, u^0, u^0, p^{-1})$ define $(u^1, p^0)$ by

$$
(3.31) \quad \left\{
\begin{array}{l}
\displaystyle \frac{\partial u^1}{\partial \tau} - v \Delta_z u^1 + v \nabla_z u^1 + u^1 \nabla_z v = \nabla_z p^0 + \nabla_x p^{-1} - \frac{\partial v}{\partial t} + v \Delta_x v + v \Delta_{xz} u^{0'} \\[1em]
\text{in } C \qquad -u^{0'} \nabla_z u^{0'} - v \nabla_x u^0 - u^0 \nabla_x v \text{ in } C \\[1em]
\displaystyle \nabla_z \cdot u^1 = - \nabla_x \cdot u^{0'} \text{ in } C \; ; \ (u^1, p^0) \text{ periodic on } C \; ; \ u^1 = 0 \text{ on } \bar{C}^2
\end{array}
\right.
$$

Note that $u^1$ always exists because $\overline{v^2} = \overline{u^0} = \overline{p^{-1}} = 0$. Note also that $u^1$ is a linear function of $u^0$ and $\nabla u^0$ and that (3.30) is automatically satisfied.

Finally, taking the mean of (3.21) on C provides us with a definition of $u^0$ :

Definition of $u^0$

For a given $(v, u^1, u^0)$ define $u^0, \bar{p}^0$ by

$$
(3.32) \quad \left\{
\begin{array}{@{}l@{}}
  \displaystyle \frac{\partial u^0}{\partial t} - v \Delta u^0 + u^0 \nabla u^0 + \frac{\partial}{\partial x_i} \overline{(v_i u^1 + u^1 v)} = \nabla_x \overline{p^0} + f - \frac{\partial}{\partial x_i} \overline{(u_i^0 \cdot u^{0'})} \\[1em]
  \nabla_x \cdot u^0 = 0 \text{ in } Q \; ; \ u^0(t=0) = u_o \; ; \ u^0|_T = u^1 .
\end{array}
\right.
$$

Remark 3.1 : Note that from (3.28) $u^0$ does not depend upon $u^0$ while from (3.31) :

$$
(3.33) \quad u^1 = \begin{array}{rcl} \overline{a} & \to & \overline{a}^0 \\ & \overline{a}^i & \to \overline{a}^i u_i^0 \\ & & \overline{a}^{ij} \to \overline{a}^{ij} \frac{\partial u_i^0}{\partial x_j} \end{array}
$$

where the a's are independent of $u^0$. Therefore the last terms in (3.32) contain
second order derivatives of $u^0$.
**Conjecture 1**

Suppose that assumptions (3.4)-(3.7) and (3.17)-(3.29) holds. Then $u^\epsilon$ converges to $u^0$ in $L^2(Q)$ and $u^\epsilon - u^1(x,t, \frac{x-u^0t}{\epsilon}, \frac{t}{\epsilon^2})$ converges to $u^0$ in $L^2(0,T;H^1(\Omega))$.

**Remark 3.2** : Equation (3.32) (with (3.33)) contains no oscillatory terms. If we could prove that the spectrum of $u^0$ is bounded when $v \to 0$ then (3.32) could be integrated on a coarse grid. Of course such a proposition is very hard to show. However if (3.31) is integrated by the method of characteristics when $v=0$ and the result is used to compute $a^0, a^i, a^{ij}$ even if $v$ is not periodic then the additional terms in (3.32) dissipates and the spectrum of $u^i$ must be finite when $v \to 0$.

**Remark 3.3** : Assumption (3.17) can be reformulated as

$$ (3.34) \quad \left| \frac{\partial v^\epsilon}{\partial t} + u^0 \nabla v^\epsilon - v \Delta v^\epsilon + v^\epsilon \nabla v^\epsilon - \nabla q^\epsilon \right| \le 0(\epsilon^{-2}) ; \quad \left| \nabla \cdot v^\epsilon \right| \le 0(\epsilon^{-1}) . $$

In this form it makes more sense than (3.17) which is sufficient to define $v$. Assumption (3.29) is due to the periodicity hypothesis.

**Remark 3.4** : If $v$ is the inertial subrange of $u$ then $u^0$, should be zero, therefore $\nabla_x q - v \Delta_x v - v \nabla_x v$ is probably small.

We end this paragraph with a situation in which the coefficients in (3.33) can be computed exactly.

**Proposition 1** : When $v \ll 1$ and $v|_{\partial z} = 0$ ($\partial z = \text{boundary of } Z$)

$$ (3.35) \quad \frac{\partial}{\partial x_i} (v_i \overline{u}^1 + \overline{u}_i v) = \frac{\partial}{\partial x_i} (\overline{v z_i} \nabla_x u^0 + u^0 \overline{\nabla_x v z_i}) - \frac{\partial}{\partial x_i} (\overline{z_i} \overline{v \nabla_x u^0}) . $$

**Remark 3.5** : The case $v|_{\partial z} = 0$ corresponds to stream lines for $v$ which do not communicate from one cell to another. Note that if $\nabla_x \cdot v = 0$ then

$$ (3.36) \quad \overline{v_j z_i} = - \overline{v_i z_j} $$

and (3.35) would further simplify into
$$ (3.37) \quad \frac{\partial}{\partial x_i} (\overline{v_i u_i}^1 + \overline{u_i v_i}^1) = \frac{\partial}{\partial x_i \partial x_j} (\overline{z_i v_j \nabla \cdot u^\circ}) $$

**Remark 3.6**: Note that (3.37) does not differ from (2.13) in that case because the stream lines of $v$ being closed $a_{ij}$ given by (2.13) are zero.

**Proof of Proposition 1 :**

Let us multiply (3.31) by $z_i e_i^j$, where $e_i^j = \delta_{jk}$, and integrate over C. Then the boundary integrals disappear and

$$ (3.38) \quad \overline{v_i u_j}^1 + \overline{u_i v_j}^1 + \overline{v_j z_i \nabla \cdot u^\circ} = \overline{v z_i \nabla \cdot u_j^\circ} + \overline{u^\circ \nabla \cdot v_j z_i} \quad i \neq j. $$

One can show as in §4 that there are no boundary layers in (3.31) when $v \to 0$. By multiplying (3.31) by $z_i e_i^i - z_j e_j^i$ one can prove that (3.38) also holds when $i=j$. (For more details see the proof of Proposition 5 in §4.)

**COMMENTS**

There are a few problems with the coupling between $u^1$ and $u^\circ$:

1.°) When $v \to 0$ it is not clear whether the terms containing $u^1$ in (3.32) dissipate energy or not.

2.°) When $v \to 0$, (3.31) becomes a hyperbolic systems but the streamlines of $v$ are not its characteristics. Therefore it is difficult to say whether (3.31) has a periodic solution for $v=0$.

3.°) It seems difficult to establish existence and unicity of $u^\circ$ given by (3.32)-(3.33) because the second order part of the operator may not be elliptic.

**CONCLUSION**

Thus an eddy viscosity-type model has been derived formally to account for the interaction between the large eddies and the small eddies. However as (3.31) and (3.32) shows the model is rather involved and possibly not well posed. Looking back into the derivation of the equations one notes that the terms which create difficulties are

i) $u^\epsilon \nabla v^\epsilon$ in (3.9) and

ii) the variable $x-u^\circ t$ in $v^\epsilon$.
Therefore in an attempt to simplify the model we drop these terms of equiva-
lently we drop $(u^e-u^o)\nabla^e$ in (3.9) and keep $x-u^o t$. We have no justification
for these simplifications except that it yields a much better system of PDE
from the mathematical view point and the above calculations can be justified
### 4. - A SIMPLER MODEL

In this section we drop the term $u^\varepsilon \nabla v^\varepsilon$ in (3.9) and we assume that $v^\varepsilon$ is not transported by the mean flow. However the results apply also to the case where $v^\varepsilon$ is transported by $u^0$ when we drop $(u^\varepsilon - u^0) \nabla v^\varepsilon$.

Therefore, given $v^\varepsilon$ and $f^\varepsilon$ we wish to study the following problem and its limit when $\varepsilon \to 0$:

$$ (4.1) \quad \left\{ \begin{array}{l} \displaystyle \frac{\partial u^\varepsilon}{\partial t} - v \Delta u^\varepsilon + (u^\varepsilon + v^\varepsilon) \nabla u^\varepsilon = \nabla p + f^\varepsilon, \quad \nabla \cdot u^\varepsilon = 0 \quad \text{in } \Omega \times ]0, T[ = Q \\[1em] \displaystyle u^\varepsilon(t=0) = 0, \quad u^\varepsilon|_{\Gamma} = 0. \end{array} \right. $$

For clarity but without loss of generality we have taken homogeneous boundary conditions.

As usual $\Omega$ is a bounded open set of $\mathbb{R}^N$ ($N=2$ or 3) with regular boundary $\Gamma$ and we will use the notations:

$$ (4.2) \quad (u,v) = \int_{\Omega} u_i v_i \, dx $$

$$ (4.3) \quad |u|_0 = (u,u)^{1/2} ; \quad |u|_1 = (\nabla u, \nabla u)^{1/2} ; $$

$$ (4.4) \quad |u|_{-1} = \sup_{v \in H_0^1(\Omega)^N} \frac{(u,v)}{|v|_1} $$

$$ (4.5) \quad H_0^1(\Omega) = \{ \phi \in L^2(\Omega) \mid \nabla \phi \in L^2(\Omega)^N ; \phi|_{\Gamma} = 0 \} $$

$$ (4.6) \quad V = \{ u \in H_0^1(\Omega)^N \mid \nabla \cdot u = 0 \} $$

$$ (4.7) \quad H = \text{closure of } V \text{ in } L^2(\Omega)^N $$

$$ (4.8) \quad L_T^2(W) = \text{space of square integrable functions of } t \in ]0, T[ $$

taking values in $W$.

In this paragraph we wish to be more rigourous mathematically.
We first show that problem (4.1) is well posed.

**Proposition 2** : If $f^\epsilon \in L_T^2(H^{-1}(\Omega)^N)$, $v^\epsilon \in L_T^2(L^4(\Omega)^N \cap H)$ and $\nabla \cdot v^\epsilon = 0$, problem (4.1) has at least one solution in the weak form

$$ (4.9) \quad \left\{ \begin{array}{ll} \left(\frac{\partial u^\epsilon}{\partial t} , w \right) + v(\nabla u^\epsilon, \nabla w) - \left( (u^\epsilon + v^\epsilon)_j \right)_j u^\epsilon, \frac{\partial w}{\partial x_j} \right) = (f^\epsilon, w) \quad \forall w \in L_T^2(V) \\ u^\epsilon \in L_T^2(V) \cap L_T^\infty(H) ; \quad u^\epsilon(t=0) = 0 . \end{array} \right. $$

Furthermore if $N=2$ or if $v^\epsilon$ and $u^\epsilon$ belong to $L^4(Q)$ then (1.9) has a unique solution.

**Proof :**

This result is a simple extension of a now classical result for the Navier-Stokes equation which can be found in LADYZHENSKAYA (1963), LIONS (1969) or TEMAM (1978) for example. The proof of the proposition follows exactly the classical proof.

Now we shall study the limit $u^o$ of $u^\epsilon$ when $\epsilon \to 0$.

### 4.1. - Case 1 : $v^\epsilon$ is bounded

**Proposition 3** : Assume that the hypothesis of proposition 1 hold. If there exists $M$ such that

$$ (4.10) \quad \sup_{(x,t) \in Q} |v^\epsilon(x,t)| \le M ; \quad \int_0^T |f^\epsilon|_{-1}^2 dt \le M \quad \forall \epsilon $$

Then there exist $v^o, u^o, f^o$ such that for a subsequence, $v^\epsilon \to v^o$ in $L^\infty(Q)$ weak star, $u^\epsilon \to u^o$ in $L_T^2(H)$ and almost everywhere in $Q$, $f^\epsilon \to f^o$ weakly in $L_T^2(H^{-1}(\Omega)^N)$ and

$$ (4.11) \quad \left\{ \begin{array}{ll} \left(\frac{\partial u^o}{\partial t}, w \right) + v(\nabla u^o, \nabla w) - \left( (u^o + v^o)_j \right)_j u^o, \frac{\partial w}{\partial x_j} \right) = (f^o, w) \quad \forall w \in L_T^2(V) \\ u^o \in L_T^2(V) ; \quad u^o(t=0) = 0. \end{array} \right. $$
**Proof :**

The existence of accumulation points $v^0, u^0, f^0$ is a direct consequence of (4.10). To obtain (4.11) it suffices to let $\epsilon \to 0$ in (4.9) and to use the a priori estimates established in the proof of proposition 1.

**Remark 4.1** : If $v^\epsilon$ depends upon $\epsilon$ by

$$ (4.12) \quad v^\epsilon(x,t) = v\left(\frac{x}{\epsilon}, \frac{t}{\epsilon}, x, t\right) \quad \forall x, t, \epsilon $$

where $v$ is a periodic function on $Y \times ]0, S[$ with respect to $y = \frac{x}{\epsilon}$ and $\tau = \frac{t}{\epsilon}$, then (BENSOUSSAN-LIONS-PAPANICOLAOU (1978))

$$ (4.13) \quad v^0(x,t) = \frac{1}{|Y|S} \int_{Y \times ]0,S[} v(y,\tau,x,t) dy d\tau \quad (|Y| = \text{volume of } y) $$

**Remark 4.2** : If $v$ has zero mean on $Y \times ]0, S[$ (which is to be expected in practice) then $v^0 = 0$ and (4.11) becomes the classical Navier-Stokes equations. This means that zero mean perturbations do not affect the large structures of the flow. We shall show that it is not the case if $v^\epsilon$ is of large amplitude.

## 4.2. Case 2 : $v^\epsilon$ is not bounded but pseudo-periodic.

Now we assume that

$$ (4.14) \quad v^\epsilon(x,t) = \frac{1}{\epsilon} v\left(\frac{x}{\epsilon}, \frac{t}{\epsilon^2}, x, t\right) $$

$$ (4.15) \quad f^\epsilon(x,t) = \frac{1}{\epsilon} f^{-1}\left(\frac{x}{\epsilon}, \frac{t}{\epsilon^2}, x, t\right) + f\left(\frac{x}{\epsilon}, \frac{t}{\epsilon^2}, x, t\right) $$

$$ (4.16) \quad \nabla_y \cdot v(y, \tau, x, t) = \nabla_x \cdot v(y, \tau, x, t) = 0 \quad \text{on } Y \times ]0, \theta[ \times Q $$

where $v, f^{-1}$, $f$ are periodic with respect to $y = \frac{x}{\epsilon}$, $\tau = \frac{t}{\epsilon^2}$, $C = Y \times ]0, S[$ and bounded in $L^\infty(C \times Q)$, i.e.

$$ Y = \bigcup_{i=1}^{N} ]0, Y_i[ \text{, } v(\dots Y_i + Y_i \dots, \tau + S, x, t) = v(\dots Y_i \dots, \tau, x, t) \quad \forall y, \tau, x, t \dots $$

As before we look for the limit of $u^\epsilon$ when $\epsilon \to 0$. Note that if $f^{-1}=0$, this limit always exists even when $v$ is not periodic. Indeed, by taking $w=u^\epsilon$ in (4.9) one finds that
$$
(4.17) \quad \frac{1}{2} \frac{\partial (u^\epsilon)^2}{\partial t} + v |\nabla u^\epsilon|^2 \le |f|_{L^\infty} |u^\epsilon|_{L^2(Q)}^2
$$

which implies the existence of a weak limit $u^o$ of $u^\epsilon$.

To find $u^o$ we look for $(u^\epsilon, p^\epsilon)$ in the form

$$
(4.18) \quad u^{\epsilon}(x,t) = u^{o}(x,t) + u^{1}\left(\frac{x}{\epsilon}, \frac{t}{\epsilon}, x, t\right) + \epsilon^{2}u^{2}\left(\frac{x}{\epsilon}, \frac{t}{\epsilon}, x, t\right) + \phi^{\epsilon}(x,t)
$$

$$
(4.19) \quad p^{\epsilon}(x,t) = p^{o}\left(\frac{x}{\epsilon}, \frac{t}{\epsilon^2}, x, t\right) + \epsilon p^{1}\left(\frac{x}{\epsilon}, \frac{t}{\epsilon^2}, x, t\right) + q^{\epsilon}(x,t)
$$

where $u^1, u^2, p^o, p^1$ are *periodic* with respect to $\frac{x}{\epsilon}, \frac{t}{\epsilon^2}$ on $C$.

Notice that

$$
(4.20) \quad \frac{\partial}{\partial x_i} \left[ u^1 \left( \frac{x}{\epsilon}, \frac{t}{\epsilon^2}, x, t \right) \right] = \frac{1}{\epsilon} \frac{\partial u^1}{\partial y_i} \Big|_{y=\frac{x}{\epsilon}} + \frac{\partial u^1}{\partial x_i} \Big|_{y=\frac{x}{\epsilon}} .
$$

Therefore when $u^\epsilon, p^\epsilon$ are replaced by (4.18) (4.19) in (4.1), the rule is to treat $y, \tau, x, t$ as independent variables provided that $\frac{\partial}{\partial \tau}$ is replaced by $\epsilon^{-2} \frac{\partial}{\partial \tau} + \frac{\partial}{\partial \tau}$, $\nabla$ by $\epsilon^{-1} \nabla_y + \nabla_x$ and $\Delta$ by $\epsilon^{-2} \Delta_y + 2\epsilon^{-1} \Delta_{xy} + \Delta_x$ where

$$
(4.21) \quad \Delta_{xy} = \sum_{i=1}^{N} \frac{\partial^2}{\partial x_i \partial y_i} .
$$

Therefore (1.1) becomes

$$
\begin{equation} \tag{4.22}
\begin{aligned}
& \left\{ \left[ \epsilon^{-2} \left( \frac{\partial}{\partial \tau} - v \Delta_y \right) - 2v \epsilon^{-1} \Delta_{xy} + \frac{\partial}{\partial t} - v \Delta_x \right] \cdot (u^o + \epsilon u^1 + \epsilon^2 u^2 + \phi^\epsilon) \right. \\
& \qquad \left. + (\epsilon^{-1} v + u^o + \epsilon u^1 + \epsilon^2 u^2 + \phi^\epsilon) (\epsilon^{-1} \nabla_y + \nabla_x) \cdot (u^o + \epsilon u^1 + \epsilon^2 u^2 + \phi^\epsilon) \right. \\
& \qquad \left. = (\epsilon^{-1} \nabla_y + \nabla_x) \cdot (p^o + \epsilon p^1 + q^\epsilon) + \epsilon^{-1} f^{-1} + f \right\} \\
& \qquad \qquad = (\epsilon^{-1} \nabla_y + \nabla_x) \cdot (u^o + \epsilon u^1 + \epsilon^2 u^2 + \phi^\epsilon) = 0.
\end{aligned}
\end{equation}
$$

$$
(4.23) \quad (\epsilon^{-1} \nabla_y + \nabla_x) \cdot (u^o + \epsilon u^1 + \epsilon^2 u^2 + \phi^\epsilon) = 0.
$$
If we write that (4.22)-(4.23) holds for all values of $(\epsilon,y,\tau,x,t)$
(instead of $(\epsilon,x,t)$) then we can identify the powers of $\epsilon$ :

terms in $\epsilon^{-1}$ in (4.22) and $\epsilon^0$ in (4.23) :

$$
(4.24) \quad \frac{\partial u^1}{\partial \tau} - v \Delta_y u^1 + v (\nabla_x u^0 + \nabla_y u^1) = \nabla_y p^0 + f^{-1}
$$

$$
(4.25) \quad \nabla_y \cdot u^1 + \nabla_x \cdot u^0 = 0 .
$$

At this point we note that (4.25) averaged over C and the periodicity of $u$
imply :

$$
(4.26) \quad \nabla_x \cdot u^0 = 0
$$

which in turn imply in (4.25)

$$
(4.27) \quad \nabla_y \cdot u^1 = 0
$$

. term in $\epsilon^0$ in (4.22) and $\epsilon$ in (4.23) :

$$
(4.28) \quad \left\{
\begin{array}{l}
\displaystyle \frac{\partial u^0}{\partial \tau} - v \Delta_x u^0 + u^0 \nabla_x u^0 + v \nabla_x u^1 + \frac{\partial u^2}{\partial \tau} - v \Delta_y u^2 + v \nabla_y u^2 - 2 v \Delta_{xy} u^1 + \\
\qquad + u^0 \nabla_y u^1 = \nabla_y p^1 + f + \nabla_x p^0
\end{array}
\right.
$$

$$
(4.29) \quad \nabla_y \cdot u^2 + \nabla_x \cdot u^1 = 0 .
$$

Taking again the mean of (4.29) on Y yields

$$
(4.30) \quad \overline{\nabla_x u^1} = 0
$$

which, from (4.24), is possible only if

$$
(4.31) \quad \nabla_x \cdot (f^{-1} \overline{y} - \overline{y} v \nabla_x u^0) = 0 .
$$
Now if we take the mean of (4.28) on C we find

$$
(4.32) \quad \frac{\partial u^o}{\partial t} - v \Delta_x u^o + u^o \nabla_x u^o + \overline{v \nabla_x u^1} = \nabla_x \overline{p}^o + \overline{f}
$$

where $\overline{\phantom{XX}}$ denotes the mean with respect to $y$ and $\tau$ on $C$. Therefore it makes
sense to define $u^1$ by (4.24), (4.27) on $C$, $u^o$ by (4.28), (4.26) on $Q$, $u^2$
by (4.28), $\nabla_y \cdot u^2 = 0$ and $\phi^\epsilon$, $q^\epsilon$ by (4.18), (4.19) then we shall show that
$\phi^\epsilon$ is small.

Definitions

Let $u^o, u^1, u^2$ together with $p^o, p^1$ be defined by :

$$
(4.33) \quad \begin{cases}
\displaystyle \frac{\partial u^o}{\partial t} - v \Delta_x u^o + u^o \nabla_x u^o + \overline{v \nabla_x u^1} = \nabla_x \overline{p}^o + \overline{f} , & \nabla_x \cdot u^o = 0 \text{ in } Q \\
u^o(t=0) = 0 , & u^o|_{\Gamma} = 0
\end{cases}
$$

where $\overline{\phantom{XX}}$ denotes the mean with respect to $y$ and $\tau$ on $C = Y \times ]0, S[$

$$
(4.34) \quad \left\{
\begin{array}{@{}l@{}}
  \displaystyle \frac{\partial u^1}{\partial \tau} - v \Delta_y u^1 + v \nabla_y u^1 = \nabla_y p^o + f^{-1} - v \nabla_x u^o, \quad \nabla_y \cdot u^1 = 0 \quad \text{in } C \\
  \phantom{\displaystyle \frac{\partial u^1}{\partial \tau} - v \Delta_y u^1 + v \nabla_y u^1} u^1 \text{ periodic on } C, \quad \overline{u^1} = 0
\end{array}
\right.
$$

$$
(4.35) \quad \left\{
\begin{array}{@{}l@{}}
  \displaystyle \frac{\partial u^2}{\partial \tau} - v \Delta_y u^2 + v \nabla_y u^2 = \nabla_y p^1 + f - \bar{f} + 2v \Delta_{xy} u^1 - u^o \nabla_y u^1 - v \nabla_x u^1 + \overline{v \nabla_x u^2} \\
  \displaystyle \qquad \nabla_y \cdot u^2 = 0 \text{ in } C \; ; \quad u^2 \text{ periodic on } C.
\end{array}
\right.
$$

Remark 4.3 : Notice that in (4.34) $u^1$ depends linearly on $\nabla_x u^o$, therefore

$$
(4.36) \quad u^1 = w^o + \sum_{i,j} w^{ij} \frac{\partial u_i^o}{\partial x_j}
$$

where $w^o, w^{ij}$ are independent of $u^o$ therefore

$$
(4.37) \quad \overline{v \nabla_x u^1} = \frac{\partial}{\partial x_k} (\overline{v_k w^o}) + \frac{\partial}{\partial x_k} (\overline{v_k w^{ij}}) \frac{\partial u_i^o}{\partial x_j} \triangleq - \nabla \cdot B - \nabla \cdot A \nabla u^o
$$

and (4.33) really means
$$
(4.38) \quad \left\{
\begin{array}{l}
\displaystyle \frac{\partial u^o}{\partial t} - \nabla \cdot (A+\nu) \nabla u^o + u^o \nabla u^o = \overline{\nabla p^o} + f + \nabla \cdot B, \quad \nabla_x \cdot u^o = 0 \quad \text{in } Q \\[1.5ex]
u^o(t=0) = 0, \quad u^o|_{\Gamma} = 0.
\end{array}
\right.
$$

**Remark 4.4** : For v, $f^{-1}$, $\nabla_x u^o$ given, $u^1$ defined by (1.34) exists and is unique.
Similarly for $u^2$, notice that in those two equations (x,t) are parameters. The
existence of $u^o$ is more difficult to establish and will be discussed later on.

**Theorem 2** : Beside (4.14)-(4.16) we assume that

$$
(4.39) \quad \overline{f^{-1}y} = 0, \quad \overline{v^y} = 0 \quad \text{(mean with respect to y only on Y)}
$$

and that $f^{-1}$ and $v$ are zero in a neighborhood of $\Gamma$ and $t=0$. If the solutions of (4.33)-(4.35) are regular ($u^o \in L^4(Q)$, $u^1, u^2 \in L^4(C \times Q)$) then there exists $K$ such that

$$
(4.40) \quad |u^{\epsilon} - u^{0}|_{0}^{2}(t) + \int_{0}^{t} |\nabla[u^{\epsilon} - (u^{0} + \epsilon u^{1})]|^{2}(s)ds \le K\epsilon.
$$

**Remark 4.5** : Conditions (4.39) are dictated by (4.31) ; $f^{-1}$, $v$ zero near $\Gamma$
and $t=0$ is required by the boundary condition of $u^{\epsilon}$ to imply $u^{1},u^{2}$ null in
a neighborhood of $\Gamma$ and $t=0$ ; the regularity of $u^{0},u^{1},u^{2}$ is for existence
of solutions and of the integrals in the proof. As in theorem 1, $w^{0}$ and $w^{ij}$
of (4.36) must be in $L^{\infty}(Q)$.

**Remark 4.6** : It is now clear that if the scaling $t/\epsilon^2$ is not introduced $u^1, u^2$
would not be solutions of parabolic systems and the development (4.18) (4.19)
could not be made.

**Remark 4.7** : The previous argument is classical in the theory of homogenization (see BENSOUSSAN-LIONS-PAPANICOLAOU (1978) for example). One can also show that $u^2$ is the corrector of order 2 to $u^o$.
**Proposition 4 (LIONS)**

Equation (4.38) has at least one solution and we have a conservation of energy in the sense :

$$
(4.41) \quad \frac{\partial}{\partial t} \frac{1}{2} |u^0|^2_0 + v |\nabla u^0|^2_0 + v |\overline{\nabla_y u^1}|^2_0 = (\bar{f}, u^0) + \int_{\Omega} \overline{f^{-1} u^1} dx
$$

Proof :

Let us multiply (4.38) by $u^0$ and integrate with respect to $x$ over $\Omega$. Then
from (4.33)

$$
(4.42) \quad \frac{\partial}{\partial t} \frac{1}{2} |u^0|^2_0 + v |\nabla u^0|^2_0 + (u^0, \overline{v \nabla_x u^1}) = (\bar{f} \cdot u^0).
$$

Now multiplying (4.34) by $u^1$ and integrating over $C$ yields

$$
(4.43) \quad \overline{v(\nabla_y u^1)^2} = \overline{f^{-1} u^1} - \overline{u^1 v \nabla_x u^0} .
$$

If (4.43) is integrated over $\Omega$ with respect to $x$ then we find that

$$
(4.44) \quad (\overline{u^0}, \overline{v \nabla_x u^1}) = - (\overline{u^1}, \overline{v \nabla_x u^0}) = v \int_{\Omega} (\overline{\nabla_y u^1})^2 dx - (\overline{f^{-1} u^1})
$$

which, together with (4.42) proves (4.41).

Therefore $-\nabla \cdot (A+\nu)\nabla$ is a strongly elliptic operator and one can show that
(4.38) has at least one solution exactly as in the case of the Navier-Stokes
equation.

Proof of theorem 2 :

We have chosen $u^0, u^1, u^2$ in such a way that (4.22) becomes

$$
\[
(4.45) \quad \left\{
\begin{aligned}
& \frac{\partial \phi^\epsilon}{\partial t} - v \Delta \phi^\epsilon + (\epsilon^{-1} v + u^0 + \epsilon u^1 + \epsilon^2 u^2 + \phi^\epsilon) \nabla \phi^\epsilon + \phi^\epsilon \nabla (u^0 + \epsilon u^1 + \epsilon^2 u^2) = \\
& = \nabla_x q^\epsilon + \epsilon (\nabla_x p^1 + 2v \Delta_{xy} u^2 + v \Delta_x (u^1 + \epsilon u^2)) - \frac{\partial u^1}{\partial t} - \epsilon \frac{\partial u^2}{\partial t} - v \nabla_x u^2 - u^0 (\nabla_y u^2 + \nabla_x (u^1 + \epsilon u^2)) \\
& \qquad - u^1 (\nabla_y (u^1 + \epsilon u^2) + \nabla_x (u^0 + \epsilon u^1 + \epsilon^2 u^2))) ; \text{ with } y = \frac{x}{\epsilon}, \tau = \frac{t}{\epsilon^2} \quad \forall (x,t) \in Q.
\end{aligned}
\right.
\]
$$

Likewise (4.23) becomes
$$
(4.46) \quad \nabla \cdot \Phi^{\epsilon} + \epsilon^2 \nabla_x \cdot u^2 = 0 \quad \text{with } y = \frac{x}{\epsilon}, \tau = \frac{t}{\epsilon^2} \quad (x,t) \in Q.
$$

Now given $\phi^\epsilon \in L_T^2(H^1(\Omega)^N)$ there exist $\eta^\epsilon \in L_T^2(H^1(\Omega)^N \cap H^2(\Omega)^N)$ and $w^\epsilon \in L_T^2(H^1(\Omega)^N)$ such that

$$
(4.47) \quad \phi^{\epsilon} = \nabla \eta^{\epsilon} + w^{\epsilon}, \quad \nabla \cdot w^{\epsilon} = 0.
$$

From (4.46) it is found that

$$
(4.48) \quad \Delta \eta^{\epsilon} = - \epsilon^2 \nabla_x \cdot u^2
$$

therefore

$$
(4.49) \quad |n^\epsilon|_1 \le K\epsilon^2 .
$$

Now we multiply (4.45) by $w^\epsilon$ and we integrate over $\Omega$

$$
(4.50) \quad \left\{
\begin{aligned}
& \frac{\partial}{\partial t} \left( \frac{1}{2} |w^\epsilon|_0^2 \right) + v |\nabla w^\epsilon|_0^2 + v (\Delta \nabla \eta^\epsilon, w^\epsilon) + \left( (\epsilon^{-1} v + u^0 + \epsilon u^1 + \epsilon^2 u^2 + \phi^\epsilon) \nabla (\nabla \eta^\epsilon), w^\epsilon \right) + \\
& \qquad \left( (\nabla \eta^\epsilon + w^\epsilon) \nabla (u^0 + \epsilon u^1 + \epsilon^2 u^2), w^\epsilon \right) = \epsilon \left( \nabla_x p^1 + 2v \Delta_{xy} u^2 + v \Delta_x (u^1 + \epsilon u^2) \right) - \frac{\partial u^1}{\partial t} - \\
& \qquad - \epsilon \frac{\partial u^2}{\partial t} - v \nabla_x u^2 - u^0 (\nabla_y u^2 + \nabla_x (u^1 + \epsilon u^2)) - u^1 (\nabla_y (u^1 + \epsilon u^2) + \nabla_x (u^0 + \epsilon u^1 + \epsilon^2 u^2)), w^\epsilon
\end{aligned}
\right.
$$

where we have used the equality : $(\nabla \eta^\epsilon, w^\epsilon) = 0$, implied by (4.47).

Recall (LADYZHENSKAYA (1963)) that

$$
(4.51) \quad |u \nabla v, w| \le \sqrt{N} |\nabla v|_0 |u|_{L^4} |w|_{L^4}
$$

$$
(4.52) \quad |u|_{L^4} \le (3\alpha |\nabla u|_0 + \frac{1}{\alpha^3} |u|_0).
$$

Then (4.50) integrated in t yields
$$
\left\{
\begin{array}{l}
\displaystyle \frac{1}{2} |w^\epsilon|^2(t) + v \int_0^t |\nabla w^\epsilon|_o^2 dt \le \int_0^t \{v\epsilon^2 |\nabla_x \cdot u^2|_o |w^\epsilon|_o + \epsilon^{-1} \sqrt{N} |\nabla \eta^\epsilon|_{L^4} \\[2ex]
\quad \times (v+\epsilon u^\circ + \epsilon u^1 + \epsilon^2 u^2 + \epsilon \nabla \eta^\epsilon)_{L^4} |\nabla w^\epsilon|_o + \sqrt{N} |\nabla \eta^\epsilon|_{L^4} |w^\epsilon|_{L^4} |\nabla w^\epsilon|_o + |\nabla \eta^\epsilon|_{L^4} \\[2ex]
\quad \times |\nabla (u^\circ + \epsilon u^1 + \epsilon^2 u^2)|_o |w^\epsilon|_{L^4} + |w^\epsilon|_{L^4} (u^\circ + \epsilon u^1 + \epsilon^2 u^2) |_{L^4} |\nabla w^\epsilon|_o + \\[2ex]
\quad + \epsilon C_1 |\nabla_x p^1| + 2v \Delta_{xy} u^2 + v \Delta x (u^1 + \epsilon u^2) - \frac{\partial u^1}{\partial t} - \epsilon \frac{\partial u^2}{\partial t} - v \nabla_x u^2 - u^\circ (\nabla_y u^2 + \nabla_x (u^1 + \epsilon u^2)) \\[2ex]
\quad - u^1 (\nabla_y (u^1 + \epsilon u^2) + \nabla_x (u^2 + \epsilon u^1 + \epsilon^2 u^2)) |_{-1} |\nabla w^\epsilon|_o \} dt.
\end{array}
\right.
$$

If $u^0, u^1, u^2$ are in $L^2(]0,T[\times]0,S[; H^1(\Omega\times Y)^N)$ then $|\frac{\partial u^1}{\partial t} + \epsilon \frac{\partial u^2}{\partial t}|_{-1}$ is finite
and we find from (4.48) and (4.52) that

$$
(4.53) \quad \frac{1}{2} |w^\epsilon|_o^2(t) + v \int_0^t |\nabla w^\epsilon|_o^2 dt \le C_2 \int_0^t \left( 3\alpha |\nabla w^\epsilon|_o + \frac{1}{\alpha^3} |\nabla w^\epsilon|_o \right) |\nabla w^\epsilon|_o dt + \epsilon C_3 \int_0^t |\nabla w^\epsilon|_o dt.
$$

Now if we choose $\alpha \le v/3C_2$ we find that it yields

$$
(4.54) \quad \int_0^t |w^\epsilon|^2 dt \le \epsilon C_4 \exp C_5 t .
$$

Therefore $w^\epsilon \to 0$ strongly in $L^2(Q)$. Now (4.54) and (4.53) imply

$$
(4.55) \quad \int_0^t |\nabla w^\epsilon|_o^2 dt \le \epsilon C_6 ;
$$

which, together with (4.49) proves the theorem.

4.3. Computation of A and B when $v \ll 1$

When $v$ is very small, $v\Delta_y u^1$ becomes negligible in (4.34) and a boundary layer
of thickness $O(v)$ may appear in the viscinity of $\partial Y$.

Form (4.34) and (4.36) we find that if $w^0$ and $w^{ij}$ remain in a bounded set of
$L^2(C)$ when $v \to 0$, then the condition $\nabla_y \cdot u^1 = 0$ implies that the limit of (4.34)
is
$$
(4.56) \quad \left\{
\begin{array}{l}
\displaystyle \frac{\partial u^1}{\partial \tau} + v \nabla_y u^1 = \nabla_y p^0 + f^{-1} - v \nabla_x u^0, \quad \nabla_y \cdot u^1 = 0 \text{ in } C \\[2ex]
\text{u}^1 \cdot n \text{ periodic on } C, \quad \bar{u}^1 = 0.
\end{array}
\right.
$$

Let us multiply (4.56) by w(y,τ) and integrate by parts on C :

$$
(4.57) \quad \left\{
\begin{aligned}
& \int_C \left( \frac{\partial u^1}{\partial \tau} w - v_i u^1 \frac{\partial w}{\partial y_i} \right) dy \, d\tau + \int_{\partial Y \times [0, \theta]} v \cdot n u^1 \cdot w \, d\Gamma \, d\tau = \\
& = \int_C (\nabla_y p^0 + f^{-1} - v \nabla_x u^0) w \, dy \, d\tau, \quad \forall w
\end{aligned}
\right.
$$

this equility implies that

$$
(4.58) \quad \int_{\partial Y \times [0, \theta]} v \cdot n \, u^1 \cdot w \, d\Gamma \, d\tau = 0 \quad \forall w
$$

which in turn implies $u^1$'s periodic or $v.n = 0$. This indicates that if (4.56)
has a solution then $u^1$ has no boundary layer (but there may be a boundary
layer for $\nabla_y u^1$ which may no longer be periodic when $v \to 0$).

If $v|_{\partial Y} = f^{-1}|_{\partial Y} = 0$ then (4.56) implies that $\nabla_y p^0|_{\partial Y} = 0$, $p^0|_{\partial Y} = 0$;
therefore taking $w = y_i e^j$, $(e_k^j = \delta_{jk})$ i≠j yields

$$
(4.59) \quad \overline{v_i u_j^1} = - \overline{f_j^{-1} y_i} + \overline{y_i v \nabla_x u_j^0} .
$$

Similarly by taking $w = y_i e^i - y_j e^j$, $i \neq j$ we find that

$$
(4.60) \quad -\overline{v_i u_i^1} + \overline{v_j u_j^1} = \overline{f_i^{-1} y_i} - \overline{f_j^{-1} y_j} - \overline{y_i v \nabla_x u_i^0} + \overline{y_j v \nabla_x u_j^0}
$$

therefore

$$
(4.61) \quad \overline{v \nabla_x u^1} = - \frac{\partial}{\partial x_i} \left( \overline{f_j^{-1} y_i} - \overline{y_i v \nabla_x u_j^0} \right) - \frac{\partial}{\partial x_i} \left( \overline{v_N u_N^1} - \overline{f_N^{-1} y_N} + \overline{y_N v \nabla_x u_N^0} \right) .
$$

However $\overline{y_i v_j} = - \overline{y_j v_i}$, therefore :
**Proposition 4** : If $v|_{\partial Y} = f^{-1}|_{\partial Y} = 0$ and $v$ is small then $u^0$ of theorem 1 is governed by

$$
(4.62) \quad \left\{
\begin{array}{l}
\displaystyle \frac{\partial u^0}{\partial t} - v \Delta u^0 + u^0 \nabla u^0 = \nabla_x p^0 + f + \bar{y} \nabla_x f^{-1} , \quad \nabla_x \cdot u^0 = 0 \quad \text{in } Q \\[2ex]
\displaystyle u^0(t=0) = 0 , \quad u^0|_{\Gamma} = 0 .
\end{array}
\right.
$$

If $v|_{\partial Y} \neq 0$ (4.56) may not have a solution and (4.34) may blow up. For example if

$$
(4.63) \quad v(x,t,y,\tau) = (\sin y_2, 0)
$$

then (4.34) becomes

$$
(4.64) \quad -v\Delta u^1 + \sin y_2 \frac{\partial u^1}{\partial y_1} = \sin y_2 \frac{\partial u^0}{\partial x_1} + \nabla_y p^0 ; \quad \nabla_y \cdot u^0 = 0
$$

which has a solution which is unbounded when $v \to 0$ :

$$
(4.65) \quad u^1 = \left( \frac{\sin y_2}{v} \frac{\partial u^0}{\partial x_i}, 0 \right)
$$
## 5. - CONCLUSION

So we have shown that under specific assumptions, to study the mean flow $u^o$ and the mean temperature $\theta^o$ one could solve:

$$ (5.1) \quad \left\{ \begin{array}{l} \displaystyle \frac{\partial u^o}{\partial t} - v \Delta u^o + u^o \nabla u^o + \nabla \cdot (A \nabla u^o) = \nabla p^o + f, \quad \nabla \cdot u^o = 0 \\[1.5ex] \displaystyle u^o(t=0) = u_o, \quad u^o|_{\Gamma} = u_1. \end{array} \right. $$

$$ (5.2) \quad \left\{ \begin{array}{l} \displaystyle \frac{\partial \theta^o}{\partial \tau} - k \Delta \theta^o + u^o \nabla \theta^o + \nabla \cdot (A' \nabla u^o) = e \\[1.5ex] \displaystyle \theta^o(t=0) = \theta_o, \quad \theta^o|_{\Gamma} = \theta_1. \end{array} \right. $$

where A and A' are symmetric tensors which are computed by integrating a PDE on a subgrid cell; A and A' actually differ only by a pressure field, entering in the calculation. An additional external field g may also be necessary in (5.1) as well as a term of the type $\nabla \cdot (Bu^o)$, if the transport of the turbulence by the mean flow is not to be neglected.

In the case of (5.2) there are evidences that A' can be calculated by averaging the second order moments of the flow field, corresponding to the homogeneous isotropic turbulence, integrated over the stream lines.

Therefore it would be wise to verify the feasibility of the present model numerically. Numerical tests are in progress, but they are not so simple since turbulence is essentially a 3-dimensional process. It seems also that some experimental data will be needed because temperature propagation in turbulence was not studied from this point of view. As a final comment, we may point out that v can always be rewritten as $\frac{1}{\epsilon} (\epsilon v^\epsilon)$ and that A can be computed from $\epsilon v^\epsilon$. Then this means that it would be wise to leave out a free parameter, $\lambda^\epsilon$, dependent of $\epsilon$ and to be adjusted to fit the experiments and replace in (5.1)-(5.2) A and A' by $\lambda^\epsilon A, \lambda^\epsilon A'$.
## ACKNOWLEDGEMENTS

We wish to thank Professors J.L. LIONS, J. MATHIEU and L.F. TSEN and Doctors B. MANTEL, J.P. QUADRAT, J.M. THOMAS, M. VIOT for their participations and valuable comments.
## REFERENCES

BENSOUSSAN A., LIONS J.L., PAPANICOLAU G. (1978)

Asymptotic methods for media with periodic structures, North-Holland

COMTE-BELLOT G., CORSIN S. (1970)

Simple criterion time correlation of full and narrow band velocity signals
in grid-generated isotropic turbulence, J. Fluid Mech., 42, p. 453.

DEARDORFF J.W. (1973)

The use of subgrid transport equation in a 3-D model of atmospheric turbulence,
J. Fluid Eng., Sept.

KOLMOGOROF A.N., (1941)

Dissipation of energy in locally isotropic turbulence, C.R. Ac. Sc. URSS, 32-16.

LADYZHENSKAYA O. (1963)

The Mathematical theory of viscous incompressible flows. Gordon and Breach.

LANDAU L., LIFCHITZ E. (1963)

Mécanique des fluides, Ed. de Moscou.

LIONS J.L., (1969)

Quelques méthodes de résolution des problèmes aux limites non linéaires.

Dunod, Paris.

MATHIEU J. (1977)

Reflexions sur quelques approches de la turbulence, Belgrade.

MEYER R.E., O'BRIEN E.E., SCOTT L.R. (1977)

Derivation of the PDE for a stochastic non-linear advection equation.

Brookhaven National Lab. Report, AMD 742.

MONIN A.S., YAGLOM A.M. (1975)

Statistical fluid mechanics of turbulence, Vol. 2

ORSZAG S.A., PATTERSON G.S. (1971)

Models and turbulence, Springer
PAPANICOLAU G.C., KOHLER W. (1974).

Asymptotic theory of mixing stochastic ordinary differential equations.

Com. pure applied Math., Vol. 27, 641-668.

REYNOLDS W.C. (1975).

Computation of turbulent flows. American Review of Fluid Mech., Vol. 8.

SCHUMANN V. (1975)

Subgrid scale model for finite difference simulations of turbulent flows in
plane channel and annuli. J. of Comp. Phys. 18, 376-404.

STRATONOVICH R.L. (1963)

Topics in the theory of random noise. Vol. 1, Gordon and Breach, New-York.

SMAGORINSKY J.S., Monthly weather Review, 91, 99-164.