Dmitry Sokolov
email: dmitry.sokolov@loria.fr

Christian Gentil
email: christian.gentil@u-bourgogne.fr

Intuitive modeling of vapourish objects

Keywords: Iterated Function System, Condensation Sets, Morphing, Cloud Modelling

Attempts to model gases in computer graphics started in the late 1970s. Since that time, there have been many approaches developed. In this paper we present a non-physical method allowing to create vapourish objects like clouds or smoky characters. The idea is to create few sketches describing the rough shape of the final vapourish object. These sketches will be used as condensation sets of Iterated Function Systems, providing intuitive control over the object. The advantages of the new method are: simplicity, good control of resulting shapes and ease of eventual object animation.

Introduction

Modeling natural phenomena such as clouds, steam, fire, smoke is a challenging task since the dawn of computer aided design systems (CAD). Standard CAD modeling approaches were developed for smooth regular shapes and are not well suited to effectively capture the intricate features of vapourish phenom-exist. We can classify it into two categories: techniques for modeling the geometry of gases and techniques for rendering scenes with gases and atmospheric effects.

The first family uses physics-based approaches, simulating the natural process that gives birth to desired shape. These methods are generally represented by dynamic systems, i.e. systems of equations that describe the evolution of a shape. As an example of such methods one can point at works of Harris et al. [START_REF] Harris | Simulation of cloud dynamics on graphics hardware[END_REF] and Miyazaki et al. [START_REF] Miyazaki | A method for modeling clouds based on atmospheric fluid dynamics[END_REF]. An animation can also be created by adjusting the parameters of the simulation system. Whereas simulations give excellent, very realistic results, they are computationally expensive and difficult to use for many animators and modelers. Setting correct physics parameters is a timeconsuming task.

The second family of methods is based on (non-physical) 3D procedural textures. A good overview can be found in a book by David Ebert [START_REF] Ebert | Texturing & Modeling: A Procedural Approach, Third Edition (The Morgan Kaufmann Series in Computer Graphics)[END_REF]. In these approaches the artist is not limited by the laws of physics, however, it can be harder to create convincing images. Textures are generated by pseudo-random noises, defining a "perturbed" 3D shape. Generally a two-level model is used: macrostructure plus microstructure. The idea is to construct a "support shape" (macrostructure) which in large represents the desired form. In the majority of cases "support shapes" are composed of spheres. Procedural modification then alters the density distribution to create the detailed wisps. That is, procedural noise functions (e.g. Perlin's noise, fractal noise, turbulence) create the microstructure to make initial "support shape" irregular. Depending on the noise function used, a number of parameters allow to adjust perturbation characteristics. Furthermore, if necessary, a specific treatment can be applied in order to simulate the nature of desired material (cloud, smoke, fire etc).

There are two ways to create animations with this family of methods. Either the artist deforms the support shape, or he modifies the microstructure. In this case the artist has to modify parameters of a pseudo-random process. The control is very delicate, notably, in order to keep the smoothness of the animation, the transformations are to be continuous. A possible solution is to generate a noise with an additional dimension corresponding to the time. For example, we can use 4D noise functions to animate a 3D cloud.

Numerous problems arise with this approach. The method perfectly fits when it is necessary to create a (static) cloud, however, it is very difficult to create a vapourish character and, especially, to animate it. The primitives become visible during the deforming the "support shape". Then, it is almost impossible to control subtle wisps on the object. For example, the method does not allow to change directly local features of the microstructure, i.e. if an artist dislikes a wisp, he can not erase it without affecting other parts of the object.

Pseudo-random nature of the microstructure makes the control very delicate and non intuitive. Artists must be familiar with such mathematical terms as "frequency", "wavelet", "Gaussian kernel" etc. Even mathematicians can fail to construct a 4D microstructure with desired characteristics.

In order to enhance the control, some authors combine the approaches [START_REF] Schpok | Modeling and animating gases with simulation features[END_REF].

Here the user disposes an interface allowing to define a "support shape", which is deformed by a simulation process.

We have developed an alternative way to create vapourish objects, giving a preference to shape control over realistic rendering. We aim to propose a very intuitive tool, a tool which can be used by any artist without special training.

The goal is to have a tool allowing to create recognizable 3D textures of gaseous objects, resembling a given draft. Figure 1 shows a smoky ship created with the method we propose. Note that the galleon is fully three-dimensional, so it can be rendered from any point of view.

To fit the constraints, we chose a model based on the fractal geometry, namely, the Iterated Function System model (IFS). The model is introduced by Hutchinson and popularized by Barnsley. In this article we generalize the notion of classical IFS, enriching it with condensation functions. This generalization inherits all the properties, established by the IFS model such as continuous dependence on parameters. As we show in the following sections, it guarantees perfect control on shapes and animations, keeping simplicity and intuition over the modeling process. The paper is organized as follows: section 2 provides the necessary background on iterated function systems. Then, sections 3 and 3.4 describe our generalization of the standard approach. Examples may be found in section 4.

Finally, section 5 concludes the paper and provides suggestions for further research.

Background

General definitions

An iterated function system or IFS consists of a finite set of transformations f i : X → X for i = 1, 2, . . . , N , where N ≥ 1 is an integer and (X, d) is a complete metric space. Let us denote it by

F = {X; f 1 , f 2 , . . . f N }. A transformation f i : X → X is strictly contractive if and only if there exists a number 0 ≤ s i < 1 such that d(f i (x), f i (y)) ≤ s i d(x, y) for all x, y ∈ X.
The number s i is called a contractivity factor for f i . If an IFS consists of strictly contractive transformations, it is called a hyperbolic IFS and the number s = max{s 1 , s 2 . . . s N } is called a contractivity factor for the IFS.

Let (X, d) be a complete metric space. Then H(X) denotes the space whose points are non-empty compact subsets of X. Given a hyperbolic IFS, let us define a contractive (with respect to the Hausdorff metric d H ) transformation F : H(X) → H(X) as follows:

F(B) = f 1 (B) ∪ f 2 (B) ∪ • • • ∪ f N (B) f or all B ∈ H(X).
Note that we use F to denote either an IFS or a corresponding transformation.

The meaning should be clear from context. In the literature the transformation F is also called the Hutchinson operator.

It is well-known [START_REF] Hutchinson | Fractals and self-similarity[END_REF] that there exists a unique fixed point A F ∈ H(X) such that A F = F(A F ). Moreover, the contraction mapping theorem [START_REF] Barnsley | Fractals everywhere[END_REF] states that for any non-empty compact subset B 0 ∈ H(X) the sequence, recurrently defined as B k+1 = F(B k ), converges to A F with respect to the Hausdorff metric as k → ∞. The resulting set does not depend on the choice of B 0 , and for this reason A F is called the set attractor of the (hyperbolic) IFS.

Deterministic algorithm

It is possible to approximate attractors in a straightforward manner by means of the contraction mapping theorem. This process is called the deterministic algorithm. Let us consider an example, where the IFS consists of three affine transformations on C (figure 2 shows an illustration):

f 1 = 1 2 z, f 2 = 1 2 z + 1 2 , f 3 = 1 2 z + 1 2 i.
The first step is to choose an arbitrary non-empty compact subset of C. Then the iteration process consists of applying each transformation f 1 , f 2 and f 3 to the current set, and then uniting the resulting sets. Then f 1 , f 2 and f 3 are to be applied to the union and so on. The limit of this process gives the attractor of the IFS, as shown above.

Random iteration algorithm

Another approach to approximate IFS attractors is the random iteration algorithm [START_REF] Barnsley | Iterated function systems and the global construction of fractals[END_REF]. An IFS with probabilities consists of an IFS together with a distribution of probabilities p 1 , p 2 , . . . , p N , positive real numbers such that p 1 +

p 2 + • • • + p N = 1.
An IFS with probabilities may be denoted as

F = {X; f 1 , f 2 , . . . , f N ; p 1 , p 2 , . . . , p N }.
The probability p i is associated with the transformation function f i for each i ∈ {1, 2, . . . , N }. If one uses the plane C for the iteration space X, then the random iteration algorithm may be expressed in pseudo-code as follows:

1 Input: An IFS with probabilities F 2 Output: An approximation (a digital image) of the set attractor A F 3 select a random point z 0 = (x 0 , y 0 ) of the plane C 4 iterate { 5 choose a transformation f i according to probabilities {p 1 , p 2 , . . . , p N }.

6 z k+1 = f i (z k ) 7 if (k > 100) draw pixel with coordinates z k+1 8 }
This schema works because the IFS is contractive: if z k is close to the attractor A F , then f i (z k ) is yet closer. Though the algorithm starts out with a random point, the distance between the solution set A F and the "dancing" point decreases exponentially. In this example we have taken 100 iterations before drawing any point. In general, the distance between A F and z 101 is inferior to the size of the pixel. Moreover, the Elton's Ergodic Theorem [6, p. 370] states that, with probability 1, the set z k is dense in A F . No fixed number of iterations is given for the algorithm. Because of the stochastic nature of the algorithm, more iterations lead to a better approximation of the exact solution.

If there are no predefined probabilities, one can choose the transformations with equal frequency in line 5 of the random iteration algorithm.

An example of an IFS with probabilities is {C; f i ; p i ; i = 1, 2, 3, 4}, where p i = 0.25 and

f 1 (z) = 1 2 z, f 2 (z) = 1 2 z + 1 2 f 3 (z) = 1 2 z + 1 2 i, f 4 (z) = 1 2 z + 1 2 + 1 2 i.
Then the random iteration algorithm proceeds as follows: an initial point, z 0 ∈ C, is chosen. One of the transformations is selected "at random" from the set {f 1 , f 2 , f 3 , f 4 }. The selected transformation is applied to z 0 to produce a new point z 1 ∈ C. Again a transformation is selected, in the same manner, independently of the previous choice, and applied to z 1 to produce a new point z 2 . The process is repeated a number of times, resulting in a finite sequence of points {z k : k = 1, 2, . . . , number of iterations}. Starting from some step, the queue of the sequence lies in the unit square with corners at (0, 0), (1, 0), [START_REF] Harris | Simulation of cloud dynamics on graphics hardware[END_REF][START_REF] Harris | Simulation of cloud dynamics on graphics hardware[END_REF] and (0, 1). If the number of iterations is sufficiently large, a picture of the filled unit square will be the result.

In this example the "dancing" point z i visits all points of the unit square The images suggest a wonderful idea. They suggest that associated with an IFS with probabilities there is a unique "density" on the attractor of the IFS.

Barnsley [8, pp. 122-129] has developed the idea of the densities using the measure theory formalism. Measures can be used to describe intricate distribution of "mass" on metric spaces. Starting from an IFS {X; f i ; p i ; i = 1, 2, . . . , N }, Barnsley introduces another IFS F = {P(X); f i ; p i ; i = 1, 2, . . . , N }, where (P(X), d P ) is the space of normalized Borel measures on X with the Monge-Kantorovich metric on it. He shows how a transformation f being contractive on X inherits this property on P(X). Therefore, the induced transformation F : P(X) → P(X) has a unique fixed point called a fractal measure or a measure attractor.

Back to our examples, both IFS have the same set attractor (the unit square), but different measure attractors.

Theoretical framework

In spite of the fact that IFS theory deals with arbitrary complete metric spaces, almost all implementations construct fractals in a 2D plane with affine contractions. Very few approaches use three dimensions or (and) non-linear transformations. For example, one can point at works by Scott Draves [START_REF] Draves | The electric sheep screen-saver: A case study in aesthetic evolution[END_REF][START_REF] Draves | The electric sheep and their dreams in high fidelity[END_REF] in the domain of rendering and special effects.

Here we propose to consider an IFS F = {H(X); f 1 , f 2 . . . , f n }, where all operators f i act on H(X), the space whose points are non-empty compact subsets of X. It is convenient to visualize the space X as a white sheet of paper, and the space H(X) as the set of all possible black ink drawings on the sheet. In order to simplify the presentation, we talk of two-dimensional images, however all the techniques may be applied to arbitrary spaces. Thus, H(X) is the set of all possible drawings on X.

In our case, induced Hutchinson operator F maps H(H(X)) onto H(H(X)).

Recall that H(X) is the space of all possible black ink drawings; then H(H(X))

(non-empty compact subsets of H(X)) is the space of all possible collections 1 of drawings. One may visualize the fixed point A F as a "book" of drawings.

IFS with condensation

A very particular case of IFS that deals with compact subsets is the IFS with (constant) condensation, introduced by Barnsley in 1988: let (X, d) be a metric space and let C ∈ H(X). Then a constant transformation f 0 : H(X) → H(X) defined by f 0 (B) = C for all B ∈ H(X) is called a condensation transformation and C is called the associated condensation set. A condensation transformation f 0 : H(X) → H(X) is a contraction mapping on the metric space (H(X), d H )

with contractivity factor equal to zero. Therefore, it possesses a unique fixed point, namely the condensation set.

Then Barnsley defines a new IFS: let {X; f 1 , f 2 , . . . , f n } be a hyperbolic IFS with contractivity factor 0 ≤ s < 1. Let f 0 : H(X) → H(X) be a condensation

1 Uncountable in general.

transformation. Then {X; f 0 , f 1 , . . . , f n } is called a hyperbolic IFS with condensation, with contractivity factor s. Note that the IFS rests in the space X. In fact, Barnsley modified the Hutchinson operator only (by adding f 0 ):

F(B) = C ∪ f 1 (B) ∪ f 2 (B) ∪ • • • ∪ f N (B) f or all B ∈ H(X).
Let us consider an example, where the IFS consists of one two-dimensional contraction f 1 (z) = 1 2 z. Clearly, the attractor consists of a single point, namely origin. But if we add a condensation transformation f 0 which transforms any compact subset to a pine tree, then the attractor becomes a series of pine trees, refer to figure 4 for an illustration.

A generalized implementation

We would like to enrich the notion of IFS with condensation. Let us consider IFS with transformations action on the space of compact subsets H(X).

Such IFS consists of contractions on (H(X), d H ). As before, a contraction f on (H(X), d H ) is a mapping with unique fixed point. This means that for any non-empty compact subset B the sequence {f k (B)}, recurrently defined by f 0 (B) = B and f k+1 (B) = f k (B) converges to the fixed point (with respect to the Hausdorff metric). In his approach Barnsley chooses a fixed point C (the pine tree in the case of figure 4). Then, clearly, the corresponding constant mapping is a contraction.

The question is: is it possible to find contraction mappings on H(X) besides obvious constant ones? Refer to figure 5 for an example. Let us suppose that there is a set C given and there is a mapping f such as any non-empty compact subset morphes to the given set C under recurrent application of the mapping.

More strictly, the Meyer's converse to Banach's theorem [START_REF] Meyers | A converse to Banach's contraction theorem[END_REF] states that if for any B ∈ H(X) the sequence {f k (B)}, recurrently defined by f 0 (B) = B and f k+1 (B) = f k (B), converges to C uniformly on a neighbourhood of C in Hausdorff metric d H , then there is a topologically equivalent metric d H on H(X)

that makes f a strict contraction (note that the random iteration algorithm works under a much weaker condition that the whole system is to be contractive on average).

In this section we build an example of such a transformation. The problem of defining contraction mappings on H(X) is closely related to the image interpolation (binary image morphing) task. Image interpolation is a general term for a set of techniques used in Computer Graphics to generate intermediate pictures between two given images. For example, Iwanowski and Serra [START_REF] Iwanowski | The morphological-affine object deformation[END_REF] propose a morphing method consisting of three steps (figure 6 illustrates the process):

1. First, for two given binary images one detects bounding boxes in order to superpose the images. Having fixed a set C ∈ H(X), the mapping

lim i→∞ f i (A) f 1 (A) = f (A) A f 2 (A) = f (f 1 (A)) f 3 (A) = f (f 2 (A)) lim j→∞ f j (B) f 1 (B) = f (B) B f 2 (B) = f (f 1 (B)) f 3 (B) = f (f 2 (B))
f B (A) = M (A, B) f ∞ B (A) = B f 2 B (A) = M (M (A, B), B) A
f C (•) := M (•, C) has a fixed point: f C (C) = M (C, C) = C.
Moreover, it is contractive with respect to the Hausdorff metric:

lim i→∞ d H (f i C (X), C) = 0 ∀X ∈ H(X).
We refer to the set C as a condensation set for a contraction mapping f C in H(X).

Let us resume: we use Iterated Function Systems of type {H(X); f 1 , f 2 . . . , f n }.

The Hutchinson operator F maps H(H(X)) onto H(H(X)). It is convenient to visualize the fixed point A F as a book of black ink drawings (recall it can be uncountable). The set attractor A F contains all fixed points of f i . Therefore,

if f i := M (•, C i )
, then it is possible to control the attractor's shape by choosing appropriate condensation sets C i .

Note that all contraction mappings in X are contraction mappings in H(X).

For example, a contractive two-dimensional affine transformation is a contraction mapping in H(X), where corresponding condensation set is a compact set, degenerated to a point. Therefore, new implementation contains the old one, inheriting all its properties such as contraction. The limit of this process gives the attractor of the IFS, as shown above.

Deterministic algorithm

Therefore, we start from a "book" with one page in it, namely bananas S.

Since there are two transformations, the number of pages in the book will double with each step. At the second step the book contains two drawings: f A (S) and f B (S). In order to clarify the notations we omit the symbol f , concatenating the indices of transformations. For example, images ABS and f A (f B (S)) are the same (we shall also call ABS an address of the image f A (f B (S))). Thus, at the second step the book contains two drawings: AS and BS. Then recurrent call of the iteration step produces four new images: AAS, BAS, ABS, BBS.

The images (sorted alphabetically) are shown in the third row of the figure 8.

Since the limit of the process does not depend on the initial choice of drawing Starting from an arbitrary non-empty compact set S (bananas, for example), the system converges to the fixed point A F . Note that the attractor does not depend on the initial set S. In few iterations only the shape of bananas disappears completely from the "book" of drawings. Note that in the case of two transformations the pages of the drawings book (being sorted by their address) represent a smooth morphing sequence between A and B. Note that in this example we have used two condensation functions, however, it is possible to use arbitrary number of transformations. Obviously, in such a case the tree will not be binary.

Random iteration algorithm

Let us take sets A and B as it is shown in figure 7. Now let us build an IFS F = {H(X); f 1 , f 2 } consisting of two contraction mappings with condensation sets A and B, built as it is shown in section 3.2. Thus, B) and X is a rectangular sheet of paper. In in this case the measure attractor µ F (rendered in figure 9) shows a uniform transition from A to B. We can approximate measure attractors using the random iteration algorithm:

f 1 = M (•, A), f 2 = M (•,
1 Input: An IFS with probabilities F = H(rectangle); f 1 , f 2 ; p 1 = 1 2 , p 2 = 1 2 , output image resolution (w, h) 2 Output: An approximation (a digital image) of the measure attractor µ Then the models were voxelized, in our example we used 256 3 voxels. Now all is ready to calculate a measure attractor of an IFS consisting of two condensation functions. Figure 14 shows a rendering result of the measure attractor. Here we imported the densities for the voxelized space to Maya fluids container to render the scene.

F A f 2 B (A) = M ( M (A, B), B) f ∞ B (A) = B f B (A) = M (A, B)
f B (A) = M (A, B) f ∞ B (A) = B f 2 B (A) = M (M (A, B), B)
The rendering stage in Maya takes about 1 minute per image (including shading, lighting, self-shadowing) and the random iteration algorithm takes less than 1 minute (at Intel Xeon CPU 1.6 GHz). Note that no hardware acceleration was used and the implementation can be greatly improved. Still, the running times are hugely inferior to the time of the modeling phase, to create such a bull mesh an experienced artist needs about one day of work. The probabilities are taken equal to 2 3 and 1 3 , respectively.

Difference between input sketches gives the effect of nebulosity. If a voxel is present in all the sketches, then it will be opaque and dense in the final render.

Otherwise it will be more transparent. In this example the probabilities were taken equal to 2 3 and 1 3 , respectively. This choice forces the final render to be closer to the first sketch. Figure 15 shows a cloudy version of the galleon we have already met in the introduction.

Since resulting images depend continuously on the parameters of transformations (input 3D models), it is easy to create animations in the usual for artists manner. Moreover, measure attractors can be imported to Maya fluids, it allows to combine our method of smoky objects creation with physics-based methods.

Conclusions

Up to our knowledge, this paper is the first attempt to use non-constant condensations for Iterated Function Systems. Besides purely aesthetic matters,

Figure 1 :

 1 Figure 1: A smoky galleon can be easily created with the new modeling tool.

Figure 2 :

 2 Figure 2: An illustration for the deterministic algorithm. Starting from any non-empty compact subset (the flower, for example), the sequence converges to the Sierpinski's Gasket.

  with equal probabilities. But what will happen if we change the probabilities? Let us consider the above IFS with the different set of probabilities {0.30, 0.20, 0.25, 0.25}. Let us also change the game rules a bit: instead of just plotting a pixel z k , we count how many times each pixel was visited by the "dancing" point z i . Once the algorithm terminates, we draw a gray-level image with darker pixels in most visited areas. Amazingly, the resulting image changes completely! Refer to figure3, where the left image shows the result for the equal probabilities, while the right one shows the result for the new set of probabilities.

Figure 3 :

 3 Figure 3: Measure attractors: different sets of probabilities lead to different distributions of "mass" on the set attractor. Left image is computed using equal probabilities {0.25, 0.25, 0.25, 0.25}, the right one is obtained with the set {0.30, 0.20, 0.25, 0.25}.

Figure 4 :

 4 Figure 4: A geometric series of pine trees, the attractor of an IFS with condensation.

Figure 5 : 7 }

 57 Figure 5: An example of a contraction mapping on (H(X), d H ). The fixed point (chosen by the user) is shown black in the final image. Any non-empty compact subset morphes to the given set by recurrent application of the mapping f .

Figure 7 :

 7 Figure 7: An example of morphing sequence. Left and right medians are shown by dark and light gray colours, respectively.

Figure 8

 8 Figure 8 shows an example of Iterated Function System made of two transformations: {H(X); f A , f B }. The condensation sets A and B are shown in the upper corners of the image. As described in the section 2.2, the first step is to choose an arbitrary non-empty compact subset. Here we have chosen a bananas image S. Then the iteration process consists of applying each transformation f A and f B to the current set, and then taking the union of the resulting sets.

SFigure 8 :

 8 Figure 8: The deterministic algorithm for an IFS made of two transformations with condensation sets A and B. Starting from an arbitrary non-empty compact set S (bananas, for

Figure 11 :

 11 Figure 11: The median set is calculated by exactly the same routine as in the previous example (refer to figure 7), but finally it is rotated by the angle, proportional to the distance between sets A and B. The measure attractor is rendered at the bottom of the figure.

A

Figure 12 :

 12 Figure 12: This example is given for another degree of rotation. The measure attractor is rendered at the bottom of the figure.

  tated median set M (A, B), however other intermediate images, e.g. M ( M (A, B), B), are completely different. Renderings of corresponding measure attractors are given at the bottom of the figures. Resulting images depend continuously on the parameters of transformations, e.g. the condensation sets and rotation degree. Therefore, by changing slightly the parameters we can obtain smooth animations.4.2. Three dimensionsUp to this moment all examples were constructed in 2D, but nothing forbids to work in 3D, let us consider an example. The task is to create a cloud creature, a bull. Here an artist created a polygonal model of the bull. Second model was created by placing randomly (polygonal) balls inside the original shape. Refer to figure 13 for the renderings.

Figure 13 :

 13 Figure 13: These quick sketches will give birth to a cloud creature.

Figure 14 :

 14 Figure 14: The measure attractor for the IFS with condensation functions, shown in figure 13.

Figure 15 :

 15 Figure 15: This galleon is created from two input models shown in the bottom row.

ena. Many scientists met the challenge and now a number of powerful methodsPreprint submitted to ElsevierDecember 16, 2015

Kluwer Academic Publishers 2000, Palo Alto (USA), 2000, pp. 81-90.

The idea is to superpose the input images in the central position (third image) and then to calculate the median image (rightmost one).

The key step is the median set calculation. The authors define the median set M (A, B) between two sets A and B as follows:

where d is the Euclidean pseudo-distance of a point to a set. In other words, the median set consists of points which lie more close to the intersection between A and B than to the complement of their union. The interpolation sequence (morphing of the first shape into the second) can be produced by the recurrent generation of new medians.

Unfortunately, this method has a lot of restrictions. First, the union A ∪ B must not be equal to X. Next, the intersection A ∩ B must not be empty. Moreover, the intersection between all connected components is to be non-empty. The main difference between the standard random iteration algorithm game and this algorithm is that at each iteration we have to store not just a point coordinates x k ∈ X, but an image of a compact subset I k ∈ H(X).

Compare the result with the bottom row of figure 8: in this case the measure attractor µ F is the "X-Ray image" of the book of drawings. We superpose all the drawings and for each point we calculate in how many drawings it appears.

Then a gray-level image is rendered using the number of occurrences.

Next example illustrates that any contraction in X is also contraction H(X), therefore, an IFS in R 2 is a particular case of IFS in H(R 2 ). Let us consider an example of an IFS with two transformations:

Bi-unit square is a square with corners at (-1, -1), (1, -1), (1, 1) and (-1, 1).

Transformation f 1 has a condensation set rendered in the left image of figure 10.

Second contraction is a non-linear transformation defined as f 2 (r, θ) = r Up to this moment we used only one type of contractions with condensations in H(X). However, there are plenty of other contraction mappings which can be constructed keeping the same condensation sets. For example, let us modify a bit the process described in section 3.


Roger Lewandowski
email: roger.lewandowski@univ-rennes1.fr

Benoît Pinier
email: benoit.pinier@inria.fr

The Kolmogorov Law of turbulence What can rigorously be proved ? Part II

Keywords: MCS Classification : 76D05, 76F65, 65M60

We recall what are the different known solutions for the incompressible Navier-Stokes Equations, in order to fix a suitable functional setting for the probabilistic frame that we use to derive turbulence models, in particular to define the mean velocity and pressure fields, the Reynolds stress and eddy viscosities. Homogeneity and isotropy are discussed within this framework and we give a mathematical proof of the famous -5/3 Kolmogorov law, which is discussed in a numerical simulation performed in a numerical box with a non trivial topography on the ground.

Introduction

We focus in this paper on the law of the -5/3, which attracted a lot of attention from the fluid mechanics community these last decades, since it is a basis for many turbulence models, such as Large Eddy Simulation models (see for instance in [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF][START_REF] Germano | Fundamentals of large eddy simulation[END_REF][START_REF] Stephen | Turbulent flows[END_REF][START_REF] Sagaut | Large eddy simulation for incompressible flows[END_REF]). Although it is usually known as the Kolmogorov law, it seems that it appears for the first time in a paper by Onsager [START_REF] Onsager | Statistical hydrodynamics[END_REF] in 1949, and not in the serie of papers published by [START_REF] Kolmogorov | The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds' Numbers[END_REF] (see in [START_REF] Tikhomirov | Selected Works of A.N. Kolmogorov[END_REF]), where the author focuses on the 2/3's law, by introducing the essential scales related to homogeneous and isotropic turbulent flows (see formula (4.13) below). In this major contribution to the field, Kolmogorov opened the way for the derivation of laws based on similarity principles such as the -5/3's law (see also in [START_REF] Chacón | Mathematical and numerical foundations of turbulence models and applications[END_REF][START_REF] Lewandowski | The Kolmogorov-Taylor Law of turbulence : what can rigorously be proved ? Handbook of applications of chaos theory[END_REF]). Roughly speaking, the -5/3's law states that in some inertial range [k 1 , k 2 ], the energy density of the flow E(k) behaves like C te k -5/3 , where k denotes the current wave number (see figure 1 This paper is divided in a theoretical part and a numerical part, in which we aim at: i) carrefully express what is the appropriate similarity assumption that must satisfy an homogeneous and istropic turbulent flow in order to derive the -5/3's law (assumptions 4.1 and 4.2 below),

ii) to theoretically derive the -5/3 law from the similarity assumption (see Theorem 4.2 below), iii) to discuss the numerical validity of such a law from a numerical simulation in a test case, using the software BENFLOW 1.0, developed at the Institute of Mathematical Research of Rennes.

Before processing items i) and ii), we discuss on different results about the Navier-Stokes equations (2.1) (NSE in what follows), that are one of the main tools in fluid mechanics, as well as the Reynolds stress (3.7) derived by taking the expectation of the NSE, once the appropriate probabilistic frame is specified. We then define the density energy E(k), which is the energy of the flow in the sphere {k = |k|} in the Fourier space. Furthermore, we introduce the concept of dimensional bases in order to properly set Assumptions 4.1 and 4.2.

The numerical simulation takes place in a computational box with a non trivial topography (see figure 3), by using the mean NSE (3.6), the k -E model (3.14), and appropriate boundary conditions supposed to model the dynamics of the atmospheric boundary layer. Atmospheric boundary layer modeling is a modern challenge because of its significance in climate change issues. We find in the literature many simulations carried out in different configurations, such as for example the case of a flat ground [START_REF] Apsley | A limited-length-scale k-epsilon model for the neutral and stably-stratified atmospheric boundary layer[END_REF][START_REF] Bou-Zeid | Large-eddy simulation of neutral atmospheric boundary layer flow over heterogeneous surfaces: Blending height and effective surface roughness[END_REF][START_REF] Fotini | Explicit filtering and reconstruction turbulence modeling for large-eddy simulation of neutral boundary layer flow[END_REF][START_REF] Porté-Agel | A scale-dependent dynamic model for large-eddy simulation: application to a neutral atmospheric boundary layer[END_REF], the case of stable or convective boundary layers [START_REF] Mason | Large-eddy simulation of the convective atmospheric boundary layer[END_REF][START_REF] Zhou | The convective boundary layer in the terra incognita[END_REF], urban simulations where building are modeled by parallelipipeds [START_REF] Michioka | Large-eddy simulation coupled to mesoscale meteorological model for gas dispersion in an urban district[END_REF], wind farms [START_REF] Portãľ-Agel | Large-eddy simulation of atmospheric boundary layer flow through wind turbines and wind farms[END_REF], realistic configurations including mountains [START_REF] Marjanovic | Investigation of model parameters for high-resolution wind energy forecasting: Case studies over simple and complex terrain[END_REF][START_REF] Zhou | Nested large-eddy simulations of the intermittently turbulent stable atmospheric boundary layer over real terrain[END_REF]. Of course, this flows is not homogeneous nor isotropic. However, the simulations shows that the curve of log 10 (E(k)) exhibits an inertial range over 4 decades, in which the regression straight line has a slope equal to -2.1424 = -5/3 (see figure 6), suggesting that the -5/3's law is not satisfied in this case.

2 About the 3D Navier Stokes equations

Framework

Let Ω ⊆ IR 3 be a C 1 bounded convex smooth domain, Γ its boundary, T ∈ IR + (enventually T = +∞), and Q = [0, T ] × Ω. The velocity of the flow is denoted by v, its pressure by p.

The incompressible Navier Stokes equation satisfied by (v, p) (NSE in the remainder) are as follows:

(2.1)

       ∂ t v + (v • ∇) v -∇ • (2νDv) + ∇p = f in Q, (i) ∇ • v = 0 in Q, (ii) v = 0 on Γ, (iii) v = v 0 at t = 0, (iv)
where v 0 is any divergence free vector fields such that v 0 • n| Γ = 0, ν > 0 denotes the kinematic viscosity, that we suppose constant for the simplicity, f is any external force (such as the gravity for example), Dv denotes the deformation tensor, ∇• the divergence operator and (v • ∇) v is the nonlinear transport term, specifically

Dv = 1 2 ∇v + ∇v t , ∇v = (∂ j v i ) 1≤ij≤3 , v = (v 1 , v 2 , v 3 ), ∂ i = ∂ ∂x i , ∇ • v = ∂ i v i , [(v • ∇) v] i = v j ∂ j v i ,
by using the Einstein summation convention. We recall that it is easily deduced from the incompressibilty condition (see [START_REF] Chacón | Mathematical and numerical foundations of turbulence models and applications[END_REF]):

(v • ∇) v = ∇ • (v ⊗ v), v ⊗ v = (v i v j ) 1≤i,j≤3 , ∇ • (2νDv) = ν∆v.
In the following, we will consider the functional spaces

(2.2) W = {v ∈ H 1 0 (Ω) 3 , ∇ • v = 0} → V = {v ∈ L 2 (Ω) 3 , v • n| Γ = 0, ∇ • v = 0},
Throughout the paper, we assume v 0 ∈ V.

Strong solutions to the NSE

Let P be the orthogonal projection L 2 (Ω) 3 → V, A and F the operators

Av = -νP ∆v, F v = P ((v • ∇) v).
By applying P to (2.1.i) in noting that P (∇p) = 0, we are led to the following initial value problem

(2.3)    dv dt = -Av + F v + P f (t), (i)
v(0) = v 0 , (ii)
where t → v(t) and t → f (t) are considered as functions valued in W and V respectively.

(2.4) ∆p = -∇ • ((v • ∇) v) + ∇ • f , derived from equation (2.1.i) by taking its divergence.
The existence of a strong solution is proved in Fujita-Kato [START_REF] Fujita | On the navier-stokes initial value problem[END_REF]. It is subject to regularity conditions regarding the initial data v 0 and the source f . The result is stated as follows.

Theorem 2.1. We assume

i) v 0 ∈ V ∩ H 1/2 (Ω) 3 , ii) f is Hölder continuous in [0, T ].
Then there exists

T = T (ν, ||v 0 || 1/2,2,Ω , ||f || C 0,α (Ω) ) such that the NSE admits a unique strong solution v = v(t). Moreover, if f = f (t, x) is Hölder continuous in Q = [0, T ] × Ω, then v(t, x), ∇v(t, x), ∆v(t, x) and ∂v(t, x)/∂t are Hölder continous in ]0, T [×Ω.
Remark 2.2. The strong solution is solution of the equation

(2.5) v(t) = e -tA v 0 - t 0 e -(t-s)A F (v(s))ds + t 0 e -(t-s)A P f (s)ds,
which is approached by the sequence (v n ) n∈I N expressed by

(2.6) v n (t) = e -tA v 0 - t 0 e -(t-s)A F (v n-1 (s))ds + t 0 e -(t-s)A P f (s)ds,
The reader is referred to [START_REF] Cannone | Harmonic analysis tools for solving the incompressible Navier-Stokes equations[END_REF][START_REF] Chemin | Wellposedness and stability results for the Navier-Stokes equations in R 3[END_REF][START_REF]Recent developments in the Navier-Stokes problem[END_REF] for more details concerning the question of strong solutions.

Turbulent solutions

Definition 2.2. We say that v is a turbulent solution of NSE

(2.1) in [0, T ] if i) v ∈ L 2 ([0, T ], W) ∩ L ∞ ([0, T ], L 2 (Ω)), ii) ∂ t v ∈ L 4/3 ([0, T ], W ) = [L 4 ([0, T ], W)] (by writing ∂ t = ∂ ∂t for the simplicity), iii) lim t→0 ||v(•, t) -v 0 (•)|| 0,2,Ω = 0, iv) ∀ w ∈ L 4 ([0, T ], W), T 0 < ∂ t v, w > dt+ T 0 Ω (v⊗v) : ∇w dxdt+ T 0 Ω ∇v : ∇w dxdt = T 0 < f , w > dt,
where for u ∈ W, F ∈ W , < F, u > denotes the duality pairing between F and u, v) v satisfies the energy inequality at each t > 0,

1 2 Ω |v(t, x)| 2 dx + ν t 0 Ω |∇v(t , x)| 2 dxdt ≤ t 0 < f , v > dt .
Remark 2.3. Once again, the pressure is not involved in this formulation. It this frame, it is recovered by the De Rham Theorem (see for instance in [START_REF] Temam | Theory and numerical analysis[END_REF]).

The existence of a turbulent solution was first proved by Leray [START_REF] Leray | Sur le mouvement d'un liquide visqueux emplissant l'espace[END_REF] in the whole space, then by Hopf [START_REF] Hopf | uber die anfangswertaufgabe fr die hydrodynamischen grundgleichungen[END_REF] in the case of a bounded domain with the no slip boundary condition, which is the case under conderation here. This existence result can be stated as follows.

Theorem 2.2. Assume that v 0 ∈ V, f ∈ L 4/3 ([0, T ], W ). Then the NSE (2.1) has a turbulent solution.

Remark 2.4. The turbulent solution is global in time, which means that it may be extended to t ∈ [0, ∞[ depending on a suitable assumption on f . However it is not known whether it is unique or not. Moreover, it is not known if the energy inequality is an equality.

The reader is also referred to [START_REF] Constantin | Navier-Stokes equations[END_REF][START_REF] Feireisl | Dynamics of viscous compressible fluids[END_REF][START_REF] Lions | Mathematical topics in fluid mechanics[END_REF][START_REF] Temam | Theory and numerical analysis[END_REF] for further results on turbulent (also weak) solutions of the NSE.

3 Mean Navier-Stokes Equations

Reynolds decomposition

Based on strong or turbulent solutions, it is known that it is possible to set a probabilistic framework in which we can decompose the velocity v and the pressure as a the sum of the statistical mean and a fluctuation, namely

(3.1) v = v + v , p = p + p .
More generally, any tensor field ψ related to the flow can be decomposed as

(3.2) ψ = ψ + ψ .
The statistical filter is linear and subject to satisfy the Reynolds rules:

∂ t ψ = ∂ t ψ, (3.3) ∇ψ = ∇ψ, (3.4)
as well as (3.5) ψ = ψ leading to ψ = 0.

We have studied in [START_REF] Chacón | Mathematical and numerical foundations of turbulence models and applications[END_REF] different examples of such filters. Historically, such a decomposition was first considered in works by Stokes [START_REF] Stokes | On the effect of the internal friction of fluids on the motion of pendulums[END_REF], Boussinesq [START_REF] Boussinesq | Essai sur la théorie des eaux courantes[END_REF], Reynolds [START_REF] Reynolds | An experimental investigation of the circumstances which determine whether the motion of water shall be direct or sinuous, and of the law of resistance in parallel channels[END_REF], Prandtl [START_REF] Prandtl | über die ausgebildeten turbulenz[END_REF], in the case of the « long time average »(see also in [START_REF] Lewandowski | Long-time turbulence model deduced from the Navier-Stokes equations[END_REF]). Later on, Taylor [START_REF] Ingram | Statistical theory of turbulence. part i-iv[END_REF], Kolmogorov [START_REF] Kolmogorov | The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds' Numbers[END_REF] and Onsager [START_REF] Onsager | Statistical hydrodynamics[END_REF] have considered such decompositions when the fields related to the flow are considered as random variables, which was one of the starting point for the development of modern probability theory.

Reynolds Stress and closure equations

We take the mean of the NSE (2.1) by using (3.3), (3.4) and (3.5). We find out the following system:

(3.6)        ∂ t v + (v • ∇) v -ν∆v + ∇p = -∇ • σ (r) + f in Q, ∇ • v = 0 in Q, v = 0 on Γ, v = v 0 at t = 0, where (3.7) σ (r) = v ⊗ v
is the Reynolds stress. The big deal in turbulence modeling is to express σ (r) in terms of averaged quantities. The most popular model is derived from the Boussinesq assumption which consists in writing:

(3.8) σ (r) = -ν t Dv + 2 3 k Id, where i) k = 1 2 tr σ (r) = 1 2 |v | 2 is the turbulent kinetic energy (TKE),
ii) ν t is an eddy viscosity.

In order to close the system, the eddy viscosity remains to be modeled. To do so, many options are avaible (see in [START_REF] Berselli | Mathematics of large eddy simulation of turbulent flows[END_REF][START_REF] Chacón | A variational finite element model for large-eddy simulations of turbulent flows[END_REF][START_REF] Chacón | Mathematical and numerical foundations of turbulence models and applications[END_REF][START_REF] Volker | Large eddy simulation of turbulent incompressible flows[END_REF][START_REF] Launder | The numerical computation of turbulent flows[END_REF][START_REF] Lewandowski | The mathematical analysis of the coupling of a turbulent kinetic energy equation to the Navier-Stokes equation with an eddy viscosity[END_REF][START_REF] Mohammadi | Analysis of the k-epsilon turbulence model[END_REF][START_REF] Sagaut | Large eddy simulation for incompressible flows[END_REF]).

One of the most popular model is the Smagorinsky's model (see for instance in [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF][START_REF] Germano | Fundamentals of large eddy simulation[END_REF][START_REF] Volker | Large eddy simulation of turbulent incompressible flows[END_REF][START_REF] Lilly | Numerical simulation and prediction of atmospheric convection[END_REF][START_REF] Stephen | Turbulent flows[END_REF][START_REF] Prandtl | Prandtl-Essentials of fluid mechanics[END_REF][START_REF] Sagaut | Large eddy simulation for incompressible flows[END_REF][START_REF] Smagorinsky | On the application of numerical methods to the solution of systems of partial differential equations arising in meteorology[END_REF][START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF]), in which (3.9)

ν t = C s δ 2 |Dv|,
where C s ≈ 0.1 or 0.2 is an universal dimensionless constant, and δ a characteristic scale, ideally the size of the smallest eddies in the flow the model is supposed to catch. This model is the fundation of the wide class of Large Eddy Simulation models. The reader will find various mathematical results concerning the Smagorinsky's model in [START_REF] Beirão | On the Ladyzhenskaya-Smagorinsky turbulence model of the Navier-Stokes equations in smooth domains. The regularity problem[END_REF][START_REF] Chacón | Mathematical and numerical foundations of turbulence models and applications[END_REF][START_REF] Volker | Large eddy simulation of turbulent incompressible flows[END_REF][START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF][START_REF] Pares | Étude mathématique et approximation numérique de quelques problèmes aux limites de la mécanique des fluides incompressibles[END_REF].

We next mention the so-called TKE model, given by (3.10)

ν t = C k √ k,
which gives accurate results for the simulation of realistic flows (see for instance [START_REF] Lewandowski | Numerical simulation of water flow around a rigid fishing net[END_REF]). In model (3.10), denotes the Prandtl mixing length, C k is a dimensionless constant that must be fixed according to experimental data. In practice, is taken to be equal to the local mesh size in a numerical simulation, and k is computed by using the closure equation (see in [START_REF] Chacón | Mathematical and numerical foundations of turbulence models and applications[END_REF][START_REF] Mohammadi | Analysis of the k-epsilon turbulence model[END_REF])

(3.11) ∂ t k + v • ∇k -∇ • (ν t ∇k) = ν t |Dv| 2 - k √ k .
The reader will find a bunch of mathematical result concerning the coupling of the TKE equation to the mean NSE in [START_REF] Brossier | Impact of the variations of the mixing length in a first order turbulent closure system[END_REF][START_REF] Bulíček | On evolutionary Navier-Stokes-Fourier type systems in three spatial dimensions[END_REF][START_REF] Chacón | Mathematical and numerical foundations of turbulence models and applications[END_REF][START_REF] Gallouët | On a turbulent system with unbounded eddy viscosities[END_REF][START_REF] Lederer | A RANS 3D model with unbounded eddy viscosities[END_REF][START_REF] Lewandowski | The mathematical analysis of the coupling of a turbulent kinetic energy equation to the Navier-Stokes equation with an eddy viscosity[END_REF].

Finally, we mention the famous k -E model that is used for the numerical simulations carried out in Section 5. In this model, E denotes the turbulent dissipation (3.12)

E = 2ν|Dv | 2 ,
and dimensional analysis leads to write

(3.13) ν t = C µ k 2 E .
The coupled system used to compute k and E is the following (see [START_REF] Chacón | Mathematical and numerical foundations of turbulence models and applications[END_REF][START_REF] Mohammadi | Analysis of the k-epsilon turbulence model[END_REF] for the derivation of these equations):

(3.14)    ∂ t k + v • ∇k -∇ • (ν t ∇k) = ν t |Dv| 2 -E . ∂ t E + v • ∇E -∇ • (ν t ∇E ) = c η k|Dv| 2 -c E E 2 k ,
where C ν = 0.09, c E = 1.92 and c η = 1.44 are dimensionless constants.

4 Law of the -5/3

The idea behind the law of the -5/3 for homogeneous and isotropic turbulence is that in the « inertial range », the energy density E = E(k) at a given point (t, x) is driven by the dissipation E . In this section, we properly define the energy density E for homogeneous and isotropic turbulent flows. We then set the frame of the dimensional bases and the similarity principle in order to rigorously derive the law of the -5/3.

Remark 4.1. For homogeneous and isotropic turbulence, one can show the identity E = 2ν|Dv | 2 = 2ν|Dv| 2 (see in [START_REF] Chacón | Mathematical and numerical foundations of turbulence models and applications[END_REF]).

Energy density of the flow

Roughly speaking, homogeneity and isotropy means that the correlations in the flows are invariant under translations and isometries (see in [START_REF] Batchelor | The theory of homogeneous turbulence[END_REF][START_REF] Chacón | Mathematical and numerical foundations of turbulence models and applications[END_REF][START_REF] Lewandowski | The Kolmogorov-Taylor Law of turbulence : what can rigorously be proved ? Handbook of applications of chaos theory[END_REF]), which we assume throughout this section, as well as the stationarity of the mean flow for simplicity. Let (4.1)

IE = 1 2 |v| 2 ,
be the total mean kinetic energy at a given point x ∈ Ω, which we not specify in what follows.

Theorem 4.1. There exists a measurable function E = E(k), defined over IR + , the integral of which over IR + is finite, and such that

(4.2) IE = ∞ 0 E(k)dk.
Proof. Let IB 2 be the two order correlation tensor expressed by:

(4.3) IB 2 = IB 2 (r) = (v i (x)v j (x + r)) 1≤i,j≤3 = (B ij (r)) 1≤i,j≤3 ,
which only depend on r by the homogeneity assumption, nor on t because of the stationarity assumption. It is worth noting that

(4.4) IE = 1 2 trIB 2 (0).
Let IB 2 denotes the Fourier transform of IB expressed by

(4.5) ∀ k ∈ IR 3 , IB 2 (k) = 1 (2π) 3 I R 3 IB 2 (r)e -i k•r dr,
We deduce from the Plancherel formula,

(4.6) ∀ r ∈ IR 3 , IB 2 (r) = 1 (2π) 3 I R 3 IB 2 (k)e i k•r dk,
which makes sense for both types of solutions to the NSE, strong or turbulent (see the section 2). It is easily checked that the isotrpoy of IB 2 in r yields the isotropy of IB 2 in k. Therefore, according to Theorem 5.1 in [START_REF] Chacón | Mathematical and numerical foundations of turbulence models and applications[END_REF] we deduce the existence of two real valued functions B d and B n of class C1 such that 1

(4.7) ∀ k ∈ IR 3 , |k| = k, IB 2 (k) = ( B d (k) -B n (k)) k ⊗ k k 2 + B n (k)I 3 .
Using formula (4.7) yields

(4.8) B ii (k) = B d (k) + 2 B n (k),
which combined with Fubini's Theorem, (4.4) and (4.6), leads to (4.9)

I R 3 B ii (k) dk = ∞ 0 |k|=k B ii (k)dσ dk = ∞ 0 4πk 2 ( B d (k) + 2 B n (k)) dk,
by noting dσ the standard measure over the sphere {|k| = k}. This proves the result, where E(k) is given by (4.10)

E(k) = k 2π 2 ( B d (k) + 2 B n (k)).
Remark 4.2. From the physical point of view, E(k) is the amount of kinetic energy in the sphere S k = {|k| = k}. As such, it is expected that E ≥ 0 in IR, and we deduce from (4.2) that E ∈ L 1 (IR + ). Unfortunately, we are not able to prove that E ≥ 0 from formula (4.10), which remains an open problem.

Dimensional bases

Only length and time are involved in this frame, since we do not consider heat transfers and the fluid is incompressible. Therefore, any field ψ related to the flow has a dimension [ψ] encoded as:

(4.11) [ψ] = (length) d (ψ) (time) dτ (ψ) ,
which we express through the couple (4.12)

D(ψ) = (d (ψ), d τ (ψ)) ∈ Q 2 . Definition 4.1. A length-time basis is a couple b = (λ, τ )
, where λ is a given constant length and τ a constant time.

Definition 4.2. Let ψ = ψ(t, x) (constant, scalar, vector, tensor...) be defined on Q = [0, T ] × Ω. Let ψ b be the dimensionless field defined by:

ψ b (t , x ) = λ -d (ψ) τ -dτ (ψ) ψ(τ t , λx ),
where

(t , x ) ∈ Q b = 0, T τ × 1 λ Ω, is dimensionless. We say that ψ b = ψ b (t , x ) is the b-dimensionless field deduced from ψ.

Kolmogorov scales

Let us consider the length-time basis b 0 = (λ 0 , τ 0 ), given by (4.13)

λ 0 = ν 3 4 E -1 4 , τ 0 = ν 1 2 E -1 2 ,
where E is the dissipation defined by (3.12) (see also Remark 4.1). The scale λ 0 is kown as the Kolmogorov scale. The important point here is that

(4.14) E b 0 = ν b 0 = 1.
Moreover, for all wave number k, and because

(4.15) D(E) = (3, -2), we get (4.16) E(k) = λ 3 0 τ -2 0 E b 0 (λ 0 k) = ν 5 4 E 1 4 E b 0 (λ 0 k),
by using (4.13). We must determine the universal profil E b 0 .

Proof of the -5/3's law

The law of the -5/3 is based on two assumptions about the flow: i) the separation of the scales (assumption 4.1 below),

ii) the similarity assumption (assumption 4.2 below).

Assumption 4.1. Let be the Prandtl mixing length. Then (4.17) λ 0 << .

Assumption 4.2.

There exists an interval

[k 1 , k 2 ] ⊂ 2π , 2π λ 0 s.t. k 1 << k 2 and on [λ 0 k 1 , λ 0 k 2 ], (4.18)
∀ b 1 = (λ 1 , τ 1 ), b 2 = (λ 2 , τ 2 ) s.t. E b 1 = E b 2 , then E b 1 = E b 2 .
Theorem 4.2. Scale separation and similarity assumptions 4.1 and 4.2 yield the existence of a constant C such that

(4.19) ∀ k ∈ [λ 0 k 1 , λ 0 k 2 ] = J r , E b 0 (k ) = C(k ) -5 3 .
Corollary 4.1. The energy spectrum satisfies the -5/3 law

(4.20) ∀ k ∈ [k 1 , k 2 ], E(k) = CE 2 3 k -5 3 ,
where C is a dimensionless constant.

Proof. Let b (α) = (α 3 λ 0 , α 2 τ 0 ).

As E b (α) = 1 = E b 0 ,
the similarity assumption yields

∀ k ∈ J r , ∀α > 0, E b (α) (k ) = E b 0 (k ),
which leads to the functional equation,

∀ k ∈ J r , ∀ α > 0, 1 α 5 E b 0 (k ) = E b 0 (α 3 k ),
whose unique solution is given by

∀ k ∈ J r , E b 0 (k ) = C(k ) -5 3 , C = k 1 λ 0 5 3 E 0 k 1 λ 0 ,
hence the result. Corollary 4.1 is a direct consequence of (4.16) combined with (4.19).

Remark 4.3. It can be shown that the law of -5/3 yields the Smagorinsky's model (3.9) (see in [START_REF] Chacón | Mathematical and numerical foundations of turbulence models and applications[END_REF]).

5 Numerical experiments

Simulation setting

The computational domain Ω is a box, the size L x × L y × L z of which is equal to (1024m, 512m, 200m) (see figure 3). The number of nodes is (256, 128, 64). The bottom of the box, plotted in figure 3, has a non trivial topography modeled by gaussian smooth domes, the height of which being equal to 50 m. We perform the simulation with ν = 2.10 -5 m 2 s -1 , which yields a Reynolds number equal to 9.10

Γ = Γ t ∪ Γ f ∪ Γ b ∪ Γ g ∪ Γ i ∪ Γ o ,
where

• Γ t is the top of the box,
• Γ f is the front face,

• Γ b is the back face,

• Γ g is the bottom of the box (the ground),

• Γ i is the inlet,

• Γ o is the outlet.

The condition on Γ i is prescribed by the Monin Obukhov similitude law [START_REF] Monin | Basic laws of turbulent mixing in the surface layer of the atmosphere[END_REF]:

(5.1) v(x, y, z, t)

| Γ i = u κ ln z + z 0 z 0 , 0, 0 t ,
where κ = 0.4 is the Von Karman constant, z denotes the distance from the ground level, the aerodynamic roughness length z 0 is equal to 0.1m, the friction velocity is expressd by:

(5.2) u = κU ref ln H ref + z 0 z 0 -
   k| Γ i = u 1/2 C -1/2 ν , E | Γ i = u 3 κ(z + z 0 )
.

On Γ b , velocity, TKE and turbulent dissipation are subject to verify the no slip and homogeneous boundary conditions, (5.4)

   v| Γg = (0, 0, 0) t , k| Γg = 0, E | Γg = 0.
On the top and lateral boundaries, we put (5.5)

   v • n = 0 on Γ t ∪ Γ b ∪ Γ f , ∇k • n = 0 on Γ t ∪ Γ b ∪ Γ f , ∇E • n = 0 on Γ t ∪ Γ b ∪ Γ f .
Finally a null gradient condition is prescribed at the outlet Γ o

(5.6)

   ∇(v • n) = 0 on Γ o , ∇k • n = 0 on Γ o , ∇E • n = 0 on Γ o .

Results

The numerical scheme we use for the simulation is based on the standard finite volume method (FVM) in space, and a Euler method for the time discretization. For the simplicity, we will not write here this technical part of the work. The reader will find comprehensive presentations of the FVM in [START_REF] Eymard | Finite volume methods[END_REF][START_REF] Joel | Computational methods for fluid dynamics[END_REF][START_REF] Jasak | Error analysis and estimation for the finite volume method with applications to fluid flows[END_REF][START_REF] Weller | A tensorial approach to computational continuum mechanics using object-oriented techniques[END_REF].

The simulation reaches a statistical equilibrium in about 180 physical seconds, which is the time at which the results are displayed. In figures 4 and 5, are plotted the values of the stramwise and spanwise components of the velocity at z = 50m, which corresponds to the dome height.

In Figure 6, we have plotted the energy spectrum of the flow at (x, y, z) = (500, 200, 50) using a log-log scale, together with a straight line whose slope is equal to -5/3 = -1, 666.... and the regression straight line of log 10 (E(k)), whose slope is about equal to -2.1424. The results call for the following comments.

i) The simulation reveals a certain reliability of the code, which suggests the convergence of the numerical method. However, the mathematical convergence of the sheme remains an open question, closely related to the question of the existence of solutions mentionned in Remark 5.1.

ii) The curve log 10 (E(k)) is an irregular curve which substancially differs from a straight line, so that we cannot conclude that numerically E(k) behaves like C te k α in some interval [k 1 , k 2 ]. Moreover, there is a gap between the slope of the regression straight line of the curve and -5/3. However, something that looks like an inertial range can be identified between k = 10 -5 m -1 and k = 10 -1 m -1 . This departure from the -5/3 law asks for the following comments and questions.

• The case under consideration yields a turbulence which is not homogeneous nor isotropic, which may explain the slope equal to -2.1424 we found.

• This simulation does not validate the Kolmogorov law or any law like E(k) ≈ C te k α . We cannot infer that such a law holds or not. Many parameters may generate the oscillations we observe in the curve log 10 (E(k)), such as any eventual numerical dissipation, a wrong choice of the constants in the k -E model which also may be not accurate, the boundary conditions we used and which may be questionable.

2 Figure 1 :

 21 Figure 1: Energy spectrum Log-Log curve

1 ,

 1 by taking U ref = 36ms -1 and H ref = 200m. The turbulent kinetic energy and turbulent dissipation are setted by(5.3)

Remark 5 . 1 .

 51 The PDE system (3.6)-(3.8)-(3.13)-(3.14) with the boundary conditions (5.1)-(5.3)-(5.4)-(5.5)-(5.6) yields a very hard mathematical problem. The existence and the uniqueness of a solution is a difficult issue, whether for global weak solutions or local time strong solutions.

Figure 2 :Figure 3 :

 23 Figure 2: Computational Box

Figure 4 :

 4 Figure 4: Streamwise direction of the flow at the z = 50m cutplane.

Figure 5 :

 5 Figure 5: Spanwise direction at the z = 50m cutplane

Figure 6 :

 6 Figure 6: Energy spectrum at the point (x, y, z) = (500, 200, 50)

  Definition 2.1. We say that v = v(t) is a strong solution to the NSE in a time interval [0, T ] if dv/dt and Av exist and are continuous in [0, T ] and (2.3.i) is satisfied there.

	Remark 2.1. In definition 2.1, the pressure is not involved. It can be reconstructed by the
	following equation

  7 . We use the mean NSE with the Boussinesq assumption, coupled to the k-E model, namely the PDE system (3.6)-(3.8)-(3.13)-(3.14). We specify in what follows the boundary conditions, by considering the following decomposition of Γ = ∂Ω:

k already denotes the TKE, and from now also the wavenumber, k = |k|. This is commonly used in turbulence modeling, although it might sometimes be confusing.