Pierre Chainais

Marc Etienne Brachet

Je Remercie Également

Sébastien Galtier

Patrice Abry

Keywords: NavierStokes, Singularities, Turbulence, Multifractal, Simulations, PIC method Détection NavierStokes, Singularités, Turbulence, Multifractal, Simulations, Méthode PIC . . . . . . . . . . . . . . . . . Wavelet transforms . . . . . . . . . . . . .

Detection and characterization of quasi singularities in turbulence

Preambule

Turbulence is a phenomenon manifesting itself in a very wide range of scale, from human size to the atmosphere and up to galaxies. At the human scale, its eects are observed in a multitude of elds including meteorology and aeronautics, but also in food industries. In some cases those eects are benecial, as it increases the mixing eciency of slowly diusive scalar. In other situations, turbulence is problematic, as it increases the energy consumption due to its high level of dissipation. All these factors makes it a subject of major interest for research, as any development can be benecial for a large range of applications.

The governing equations for uids mechanic in general and turbulence in particular are the NavierStokes equations. Despite having been derived in the 18th century, they are still a subject of active research. A major yet unsolved issue is the regularity of its solutions. This issue has physical implications. Indeed, it was proven that if singularities exist in turbulent velocity elds, they would be able to dissipate energy in a non viscous way. This contribution to the global energy budget might be an explication to the phenomenon of anomalous dissipation, or the fact that the energy dissipation rate no longer depends on viscosity for high enough Reynolds numbers.

The objective of this thesis is to build a new tool to detect potential singularities in both numerical and experimental turbulent velocity elds. This tool is based on a measure of local regularity inherited from the concept of multifractal analysis.

It is meant to complement a previously designed criterion based on local energy transfer. Indeed, this other tool is designed to analyze the "anomalous dissipation" part of the problem and does not provide a way to quantify the regularity of a given velocity eld. The conjunction of these two tools allows for an in-depth analysis of the extreme events found in turbulent elds. This analysis brings some hints about the distribution, physical mechanism and evolution of the events assimilated to quasi singularities.

In parallel with the elaboration of tools to characterize the ow regularity, we develop a new numerical simulation method based on the Clebsch decomposition. 1

Chapter1

Turbulence and singularities

In this chapter, we introduce the properties of turbulence and the NavierStokes equation and show why there are interrogations about the existence or non existence of singularities.

Generalities about turbulence

Turbulence manifests itself in a ow when the viscosity of the uid is low enough that viscous eects are dominated by inertia. In this section, we recall the governing equations and present some of their properties. We particularly focus on properties leading to the suspicion of singularities.

NavierStokes: A governing equation

The governing equations for viscous ow are the Navier-Stokes equations. Derived in the 19th century, they include two partial dierential equations, one coming from the conservation of mass and the other from the conservation of momentum.

In the incompressible case, they are expressed as follow:

∂ t u i + u j ∂ j u i = - 1 ρ ∂ i p + ν∂ j ∂ j u i + f i ,
∂ j u j = 0.

(1.1)

(1.2)
where the Einstein convention for repeated indices is used.

In these equations, hereafter referred to as INSE for Incompressible Navier Stokes Equations, u i (x, t) for i = 1, 2, 3 are the three components of the velocity eld; p(x, t) is the pressure eld; ν is the molecular viscosity; f i (x, t) corresponds to the volumic forces acting on the uid.

The rst equation is derived from the conservation of momentum (i.e. Newton's second law) applied to an elementary volume of uid. In this equation, the left hand side, corresponding to a material derivative of the velocity ( du i dt ), is the elementary volume acceleration. The right hand side corresponds to the forces applied to the volume. This includes a pressure term, the contribution of viscosity, as well as exterior volumic forces. The second equation expresses mass conservation, or in this incompressible case, the conservation of volume.

Several conditions must be fullled for the INSE to be applicable to a ow. The rst of those is the incompressibility condition. Most uids can be considered incompressible, i.e. keeping an uniform density ρ, under certain conditions. For this to be valid, the time required to regularise variations of density must be very small compared to the other characteristics time of the ow. This property is captured by the Mach number:

M = u c (1.3)
where u is a characteristic velocity of the ow and c is the sound velocity in the medium. The incompressibility condition is then valid provided M 1.

Another condition is the scale separation, related to the validity of using an elementary volume of uid when computing the conservation of mass and the conservation of momentum. This is made more apparent when the NavierStokes equations are derived using the kinetic theory of gas (see e.g. [START_REF] Lachowicz | From kinetic to Navier-Stokes-type equations[END_REF]). The general idea is that for the momentum equation (1.1) to be valid without including molecular eects, the elementary volume must be very large compared to the range of these interactions. This property is captured by the Knudsen number:

Kn = mol min,f low (1.4)
where mol is the scale of the molecular interactions and min,f low is the minimum scale of the ow. This last scale represents the size of the elementary volume of uid that needs to be considered in order to describe the ow with equation (1.1).

Typically, min,f low is the Kolmogorov scale. The scale mol is typically the mean free path for a gas or the distance between molecules for a liquid. In this thesis, we only use uids with a molecular scale of the order of 0.1 nm, while the Kolmogorov scale is always greater than 10 µm. As a result, the scale separation should always be respected.

Inertia versus viscosity

The momentum part of the INSE can be made dimensionless using a characteristic velocity U and a characteristic scale L as follows:

∂ * t u * i + u * j ∂ * j u * i = -∂ * i p * + 1 Re ∂ * j ∂ * j u * i + f * i (1.5)
where the asterisk denotes a dimensionless quantity or derivative and the Reynolds number is dened as:

Re = U L ν .

(1.6)

The Reynolds number is thus the main parameter for the INSE. It can be interpreted as a ratio of inertial eect over viscous eects. Indeed, if the Reynolds number is low, the viscous term ∂ * j ∂ * j u * i dominates the inertial terms in the INSE.

This leads to a laminar ow, with a number of properties including symmetry by time reversal. Such ows can be observed for very viscous uid (e.g. honey) or at very small scale and velocity (e.g. ow around a bacteria swimming).

Conversely, if the Reynolds number is very high, the nonlinear terms u * j ∂ * j u * i and ∂ * i p * dominate the viscous term. This type of ow is much more chaotic than laminar ows, leading to the name turbulence, originating from the Latin turba meaning "disorder" or "tumult". Such ow is commonly seen in the wakes of boats or planes for example.

In practice, the exact threshold of Reynolds number separating laminar and turbulent ow depends on the geometry of the system.

A multiscale phenomenon

One characteristic of turbulent ows is the corresponding wide range of scales. This property is observed visually in the velocity elds as one can observe vortices of very dierent size in a same ow. Examples of such velocity eld are given in Chapter 4.

Meanwhile, the multiplicity of scales is also something that can be anticipated from the INSE. This is particularly the contribution of the non linear terms, fundamental to turbulence. Indeed, let us suppose an initial condition of velocity with only one mode of wave number k. By the trigonometric formula cos(kx) 2 = cos(2kx) + 1 2 ,

(1. 7) this initial mode will cause the apparition of a mode with wave number 2k. This process is repeated, leading to cover a large range of scales. In the limit of zero viscosity, corresponding to the Euler equations, this cascade continues indenitely.

However, in the case of a nite viscosity, this process is interrupted as the viscous term becomes comparable to the non linear term. Using the dissipation rate as a basis for dimensional analysis, one can estimate the characteristic velocity at scale as u = ( ) 1/3 . This means that the corresponding viscous term (resp. nonlinear term) has an order of magnitude of ν 1/3 -5/3 (resp. 2/3 -1/3 ). Those terms are balanced at the scale = η, known as the Kolmogorov scale:

η = ν 3 1/4 . (1.8)
This "cascade" can be observed in the energy spectra corresponding to turbulent velocity elds. Such spectra are shown in Chapter 3, and a schematic illustration is presented in gure 1.1. The main feature observed is a large range of modes, named the inertial range, over which the energy follows a power law in k -5/3 . This is the Kolmogorov law, named after Andrey Kolmogorov who derived it in his famous 1941 paper [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous uid for very large Reynolds numbers[END_REF]. It can be recovered less rigorously using an argument based on dimensional analysis. Indeed, as the inertial range observed experimentally occurs above the Kolmogorov scale, the viscous eect are still dominated by inertial eects. Removing viscosity, the only pertinent physical parameter for the analysis is the energy dissipation rate (in a stationary regime, this is also the energy injection rate). As such, the energy of the mode of wave number k must scale as E(k) ∝ 2/3 k -5/3 . The inertial range is followed at small scale by the dissipative range (also referred to as the viscous range). For scales close to the Kolmogorov scale, viscosity becomes eective and energy is dissipated. On the other side of the inertial range, at large scale, is the injection range. This corresponds to the scales at which energy is injected in the ow. The exact behavior at this range is highly dependent on the method of forcing and the geometry of the ow.

The zeroth law of turbulence

The viscous term in the INSE leads to energy dissipation. The energy balance assuming smooth solutions reads:

∂ t 1 2 |u| 2 + ∇ • 1 2 |u| 2 + p u -ν∇ 1 2 |u| 2 = -ν|∇u| 2 (1.9)
where the dissipative term is on the right hand side. Thus, the average dissipative term = ν|∇u| 2 , the average being taken in space and time, appears to scale with the viscosity ν. In other world, it would be expected that the lower the viscosity of the uid (i.e. the higher the Reynolds number), the less forcing is required to keep it in motion. This intuition is however disproved by both experiment and simulations. In what is referred to as the zeroth law of turbulence or anomalous dissipation. This behavior is schematically illustrated in the gure 1.2. What is observed is that for low values of Reynolds number, the dissipation rate indeed scales proportionally to the inverse of Re, i.e. to the viscosity. However, past a certain critical Reynolds number, this scaling changes and the dissipation rate no longer depends on the Reynolds number [START_REF] Pearson | Measurements of the turbulent energy dissipation rate[END_REF].

log Re log ∝ Re -1 ∝ Re 0 Reynolds number leads to a rst argument for the singular behavior of the ow.

Indeed, if the term becomes independent of ν for high Reynolds number, this means that the gradients of velocity ∇u become larger and larger as the viscosity increases. As the corresponding scaling would be ∇u ∝ Re -1/2 , the limit Re → ∞ leads to a divergence of the velocity gradient.

From another point of view, at ν = 0 the NavierStokes equations 1.1 become the Euler equation. This equation is such that, assuming smooth solutions, the energy is conserved. As such, it appears paradoxical that the limit ν → 0 is observed to have a constant, non zero dissipation. Lars Onsager [START_REF] Onsager | Statistical hydrodynamics[END_REF] proposed a solution to this contradiction: a non smooth velocity eld solution of the Euler equations would allow for a non zero dissipation. . . . in three dimensions a mechanism for complete dissipation of all kinetic energy, even without the aid of viscosity, is available." (Onsager)

Mathematical results

In this thesis, we are interested in looking for such loss of smoothness. In particular, a major problem which has yet to be resolved since the 19th century is whether the 3D incompressible Euler and NavierStokes equations are well posed. The later in particular is the object of one of the "Millenium Problems" set by the Clay Mathematics Institute of Cambridge. This problems are a collection of dicult mathematical questions, the resolution of which would be rewarded by a prize of one million dollars.

The well posedness of an equation has three requirements in the meaning of Hadamard [START_REF] Hadamard | Sur les problèmes aux dérivées partielles et leur signication physique[END_REF]. Those are the existence of a solution, its unicity and its regularity.

For this last criterion, it appears indeed that if a solution to NavierStokes is not at least twice dierentiable, the Laplacian in the viscous term of equation (1.1) is not properly dened.

In the 2D case, the well posedness of the Euler equations [START_REF] Yudovich | Non-stationary ows of an ideal incompressible uid[END_REF] and the Navier Stokes equations [START_REF] Ladyzhenskaya | The mathematical theory of viscous incompressible ow[END_REF] have already been proven. The proof in the 3D case is much more dicult. For the 3D incompressible Euler equations, it was proven that there is no unicity of the solution, if it exists [START_REF] Shnirelman | On the nonuniqueness of weak solution of the Euler equation[END_REF]. Moreover, it was proven that if a non regular solution exists, then the vorticity is unbounded [7]. There are however no result proving the existence of solutions, nor their regularity or non regularity.

For the case of the Navier Stokes equations, Leray [START_REF] Leray | Sur le mouvement d'un liquide visqueux emplissant l'espace[END_REF] actually proved the existence of solutions in the weak sense, without any assumption about regularity. This leaves uniqueness and regularity as open questions. Some results have been obtained for this second point. First, Constantin [14] proved that if singularity exists in NavierStokes, then they take the form of an unbounded velocity |u| → ∞. This is an even stronger condition than for the Euler case. Moreover, it was also proven that singularity must be very rare events. More precisely, singularities cannot form continuous curves in space-time, i.e. they must be punctual, instantaneous events [11]. Finally, singular stationary solutions of NavierStokes have been found [START_REF] Li | Homogeneous solutions of stationary NavierStokes equations with isolated singularities on the unit sphere[END_REF]. This however does not provide a proof for the regularity of solutions starting from a regular initial condition, which would be more pertinent from a physical point of view.

Physical approach to singularities

The impact of those mathematical results on our understanding of the physical ow described by the NavierStokes equations is not obvious. In particular, the physical reality of a point like singularity or a diverging velocity is questionable.

Under these circumstances, it is important to make the distinction between the mathematical equations and the phenomena observable in real ows, as well as simulated ows.

The objective of this thesis is to approach the question of the regularity of the solutions to the NavierStokes equations from a physical point of view. The working hypothesis is that mathematical singularities do exist, i.e. that there exist solutions becoming singular in a nite time from a regular initial condition. This hypothesis does not lead to the existence of singularities in a real ow. Indeed, for the INSE to work as a model of the ow, the conditions exposed in section 1.1.1 must be fullled. However, the very concept of punctual singularity, in other words a structure at scale = 0, is contradictory with the condition of scale separation.

Furthermore, the mathematical result by Constantin [14] means that a singularity of Navier-Stokes would lead to an unbounded velocity. From the strict perspective of the NavierStokes equations, this would mean breaking the incompressibility condition, as the Mach number would locally blow up: M a 1.

In summary, the existence of singularities in NavierStokes would lead to locations of the ow where the INSE are no longer the valid model. The detection of such events, hereafter called potential singularities or quasi singularities, is the target of this thesis. In accordance to the mathematical results already presented [11], such events are expected to be rare. This means that there is no reason to question the validity of the NavierStokes equations in general to describe a turbulent ow. However, in accordance with Onsager's conjecture, the mathematical singularities could lead to a non viscous energy dissipation, enough to explain anomalous dissipation. This would mean that quasi singularities might be relevant to the global energy budget despite their rarity. The questions of energy transfers are further discussed in the next chapter.

In this thesis, we will primarily work using numerical simulations, even though all tools described are similarly applicable to experimental results. By their nature, simulations fundamentally dier from the real ow in the relationship with NavierStokes singularities. Indeed, as the simulations are based on the equations, the conditions of validity of the INSE do not matter. The dierence with the complete NavierStokes equations manifests itself in the numerical scheme used to discretize them. For example, a spatial derivative can be approximated using a nite dierence method using the following method (in 1D): df dx (x) ≈ f (x + ∆x)f (x -∆x) 2∆x .

(1.10)

These approximations fundamentally relies on the Taylor expansion of the function f . In other words, a certain regularity is required for the derivatives to be accurate. On the other hand, the spectral method uses exact derivatives of functions decomposed over Fourier modes. The approximation in the spectral method comes from the limited number of modes that can be computed, eectively cutting o uctuations at a lower scale. In other words, the accuracy of the simulation of the INSE requires a grid step smaller than the smallest scales of the ow, such that ∆x < min,f low . Typically, the discretization scale ∆x is chosen of the order of the Kolmogorov length η. This means that if the solution to the NavierStokes equations contains real mathematical singularities, those cannot be observed in the equivalent simulations. Note that more complex methods, for example using an adaptive mesh, may not have a denite smallest scale. Even for these methods, the scale = 0 of a punctual singularity is still inaccessible.

We expect that wherever a real singularity would happen in the INSE, we would detect a quasi singularity in the simulated velocity eld, resolved down to ∆x.

This justies our use of simulations for the detection of singular features. In this aspect, our approach is similar to several previous studies which exhibited emergent irregular structures in simulations. These studies are preferentially conducted using the Euler equations, like the works by Luo and Hou [START_REF] Luo | Potentially singular solutions of the 3D axisymmetric Euler equations[END_REF] or Pumir [START_REF] Pumir | Development of singular solutions to the axisymmetric euler equations[END_REF],

who observed events with a diverging vorticity. Similarly, numerous numerical studies have been performed by Kerr (e.g. [START_REF] Kerr | Evidence for a singularity of the three-dimensional, incompressible euler equations[END_REF][START_REF] Kerr | The role of singularities in turbulence[END_REF][START_REF] Kerr | The role of singularities in Euler[END_REF][START_REF] Kerr | Vortex collapse and turbulence[END_REF]) in order to study the relationship between potential singularities and vortex dynamics using criteria based on enstrophy and vorticity. Other such simulations include works by Grauer et al. [START_REF] Grauer | Finite time singularities in ideal uids with swirl[END_REF], Uhlig et al. [START_REF] Uhlig | Singularities in cascade models of the Euler equation[END_REF] and Grafke et al. [START_REF] Grafke | Numerical simulations of possible nite time singularities in the incompressible Euler equations: comparison of numerical methods[END_REF]. Some studies have been performed using the full NavierStokes equations. This is the case in particular of a work by Vassilicos et al. [START_REF] Vassilicos | Self-similar spiral ow structure in low Reynolds number isotropic and decaying turbulence[END_REF], which uses a criterion based on the geometry of the streamlines of velocity to search for singularities. This thesis uses Navier Stokes simulations and a new detection tool in order to locate and quantify quasi singularities.

We discussed that the real ow and the simulated one would no longer verify the INSE in case of a mathematical NavierStokes singularity for dierent reasons. As such, one might worry that the quasi singularities observed in these two cases would be dierent. However, we will suppose that the INSE represent a valid model down to the Kolmogorov scale even in presence of quasi singularities. Likewise, we expect the INSE to be simulated accurately down to the same scale. In consequence, the events qualied as potential singularities should be the same in experimental and simulation data. In this thesis, only simulations will be studied but some element of comparison with the experimental results by Paul Debue [15] are provided.

Additionally, the tools developed are designed to be later used on experimental data.

Chapter 2

Tracking singularities using local energy transfer and the turbulent force

As the anomalous dissipation is suspected to be related with singularities, it is natural to try to use an approach based on energy transfer when searching for them. As a result, local energy budgets have already been used as the main criterion to detect and quantify singularities, in particular in the thesis by Dennis Kuzzay [START_REF] Kuzzay | Investigations on the relevance of Onsager's conjecture in real incompressible turbulence[END_REF] and Paul Debue [15]. More precisely, the criterion is based on a scale by scale energy budget. The underlying idea is that if singularities exist and are related to the anomalous energy dissipation, then they would transfer energy to lower and lower scales without ever being dissipated by viscosity. This would go on until some hypothesis required for the formalism is broken and the energy is dissipated in a non viscous way. In contrast, in absence of singularity, the cascade of energy to lower scales would be interrupted by viscous dissipation. As such, locating regions of extreme energy transfer at a dissipative scale is an intuitive way to probe for singularities. This method is also used in this thesis. As a result, we detail the derivation of the local energy budget and some previous results in this chapter. We also describe a similar method based on the turbulent force.

2.1 Derivation using the weak formulation

Smoothed velocity eld

The formalism used was introduced by Duchon and Robert in their paper [START_REF] Duchon | Dissipation d'énergie pour des solutions faibles des équations d'Euler et Navier-Stokes incompressibles[END_REF].

The derivation of this formalism is based on an energy budget on the weak INSE. This is akin to smooth the velocity elds using a regularizing function, such that supposing their regularity is no longer required. Said another way, for all terms 14CHAPTER 2. Tracking singularities using local energy transfer and the turbulent force involving a derivative in the equation, the derivative operator is applied to the regularizing function instead of the velocity eld. The smoothed velocity eld u in 3D for a scale is expressed as a convolution of the velocity eld and the smoothing function φ :

u (x) = φ * u(x) = φ (y)u(x -y)dx. (2.1)
The regularizing function φ is dened by:

φ (x) = φ(x/ ) φ(y/ )d 3 y (2.2)
where φ is an even function, positive of norm unity. As a consequence φ is also even, positive with norm unity. As a result of the ltering, increasing the scale reduces the amount of small scales in u . On the other hand, we can verify that in the sense of distributions, lim →0 u = u.

In this thesis, for all computation related to the Duchon Robert formalism, we will use a simple Gaussian:

φ(x) = 1 (2π) 3/2 exp - x 2 2 .
(2.3) Spatial derivatives of the ltered velocity eld u can thus be interpreted as wavelet transforms of the velocity using a Gaussian wavelet.

Energy budget on the ltered velocity by Duchon and Robert

Applying the convolution with φ to the INSE (without forcing) yields:

∂ j u j =0, ∂ t u i + (u j ∂ j u i ) = -∂ i p + ν∂ j ∂ j u i . (2.4) (2.5)
By computing an energy budget on these equations, Duchon and Robert [START_REF] Duchon | Dissipation d'énergie pour des solutions faibles des équations d'Euler et Navier-Stokes incompressibles[END_REF] obtained the following relation:

1 2 ∂ t u i u i + ∂ i T i = 1 2 ν∂ j ∂ j (u i u i ) -ν∂ j u i ∂ j u i -D I (2.6)
where:

T i = 1 2 u j u j u i + p u i + pu i + 1 4 (u i u j u j ) -(u j u j ) u i (2.7)
and

D I (x) = 1 4 ∇φ (ξ) • δ ξ u(x)(δ ξ u(x)) 2 dξ. (2.8)
where δ ξ u(x) = u(x + ξ)u(x) is the velocity increment.

They then demonstrated that in the limit → 0, the equation (2.6) converges in the sense of distributions to:

∂ t u 2 i 2 + ∂ j u j u 2 i 2 + p = ν∆ u 2 i 2 -ν∂ j u i ∂ j u i + D(u) (2.9)
where D(u) = lim →0 D I . This limits holds for all smoothing function φ verifying the properties stated in the previous section.

This result adds the term D(u) to the classical energy budget. If the velocity eld is regular, a Taylor expansion leads to δ ξ u ≈ ξ • ∇u for small increments. In consequence, the scaling of the local energy transfer is D I ∝ 2 and hence we have D(u) = 0. This is however not true in the presence of strong enough singularities, which would lead to an additional energy dissipation independent of viscosity. This property is also true in the case of the Euler equations. This term is proven positive for the weak solutions of the INSE proposed by Leray.

Physical interpretation

For our objective of evaluating energy transfers in situations where the velocity eld might not be regular, the equation (2.6) is still problematic as computing the right hand side requires to apply derivative on the unsmoothed velocity eld u.

We can rearrange the terms to obtain this equation [START_REF] Dubrulle | Beyond Kolmogorov cascades[END_REF]:

1 2 ∂ t u i u i + ∂ i J i = -D ν -D I (2.10)
where

J i = T i - 1 2 ν∂ i [(u j u j ) + u j u j ] + νu j ∂ i u j (2.11) and D ν (x) = ν ∇ 2 φ (ξ) δ ξ u(x) 2 2
dξ.

(2.12)

We can interpret the dierent terms of the equation as follows:

the spatial transport, pressure term and viscous diusion are contained in the term ∂ i J i .

the viscous dissipation at scale : D ν .

16CHAPTER 2. Tracking singularities using local energy transfer and the turbulent force the energy transfer to the lower scales: D I .

Using this interpretation, we can dene a strategy on how to detect potential singularities. Indeed, if it exists, a true singularity would carry energy to the scale → 0. Similarly, a quasi singularity can be expected to transport energy to a scale smaller than the Kolmogorov scale η. As a result, the general idea of using local energy transfers to detect rare irregular structures consists in nding locations where the term D I stays high at low scale. The natural quantity to which compare D I is the viscous term D ν . As such, the method to detect quasi singularities is to identify the locations where the energy transfer D I is greater than the energy dissipation D ν at a scale ∼ η.

In order to further rene this intuition, let us suppose the velocity eld is such that velocity increments over a distance verify a scaling δ (u) ∼ h , with 0 < h < 1. This concept is formalized in the next chapter. By injecting this scaling in equations (2.8) and (2.12), we get the scalings:

D I (x) ∝ 3h-1 , D ν (x) ∝ ν 2h-2 . (2.13) (2.14)
As a result, we can obtain an estimate for the scale at which the energy transfer term D I and the energy dissipation term D ν are balanced:

η h ∝ ν 1/(h+1) . (2.15)
Those two terms being balanced means that the velocity eld starts to get regularized at scale η h . In other words, the velocity increments no longer follow the initial power law δ (u) ∼ h for < η h .

For a Hölder exponent of h = 1/3, which would correspond to the K41 scaling [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous uid for very large Reynolds numbers[END_REF] in the inertial range, this regularizing scale matches the classical Kolmogorov scale η h=1/3 = η. However, if the velocity eld is more singular, i.e. h < 1/3, this means that the regularizing scale η h gets lower than η. In the case of a true singularity, we would have η h = 0 as the velocity eld is never regularized, leading to a Hölder coecient of h = -1. The case h < 0 has not been discussed in section 3.1 and would correspond to a local blowup of velocity, which is coherent with the result by Constantin [14].

Using those scalings, our intuition of the energy budget around a potential singularities can be rened: For a singularity of "regularity" h < 1/3, the dissipation term D ν only catches up to the energy transfer term D I below the Kolmogorov scale. However, this interpretation is only based on scalings and might not be perfectly robust. In order to conrm the strength of a singularity located using this criterion, it might be necessary to actually extract the power law behavior of these terms with respect to . The computational cost of this operation makes it impractical. A more easily implemented choice is to consider any event with an extreme value of energy transfer D I as a candidate for singularity [START_REF] Kuzzay | New method for detecting singularities in experimental incompressible ows[END_REF]. This situation is one of the reason why the aim of this thesis is to develop a method to compute local Hölder exponent that could be used as a direct criterion for the detection of quasi singularities.

Turbulent force

We can also study the ltered INSE (2.5) for the velocity. Following a process much similar to the previous one, the following equation was derived by Eyink [START_REF] Eyink | Turbulent cascade of circulations[END_REF]:

∂ t u i + u j ∂ j u i = -∂ i p + ν∂ j ∂ j u i + F i (2.16)
where F can be identied as a force perceived at scale caused by the uctuations of the ow at smaller scales. This force will hereafter be referred to as the turbulent force and is computed as follow:

F i (x) = 1 ∂ j φ (ξ)δ ξ u i (x)δ ξ u j (x)dξ -∂ j φ (ξ)δ ξ u i (x)dξ φ (ξ )δ ξ u j (x)dξ .
(2.17)

One might interpret this turbulent force from the point of view of the Large Eddy Simulations (LES). In the context of LES, the turbulent force would correspond to the gradient of the Reynolds tensor. This term is then modeled, in order to avoid resolving the small scales required for its computation.

Following the same method as for the local energy transfer, one can evaluate the scaling of the turbulent force dened in equation (2.17) as:

F ∝ 2h-1 . (2.18)
This means that the norm of the turbulent force diverges when → 0 whenever the local h exponent veries h < 1/2. In other words, a criterion based on the anomalous turbulent force would be more sensitive than one based on the local energy transfer. 18CHAPTER 2. Tracking singularities using local energy transfer and the turbulent force

Implementation

The explicit computation of the integrals in the energy transfer from equations (2.8), the viscous dissipation from equation (2.12), and the turbulent force from equation (2.17) would be too computationally intensive. It is however possible to rearrange the terms in a more favorable way.

Using the Einstein summation convention, the equation (2.8) can be rewritten as:

D I (x) = 1 4 (∂ i φ )(ξ) • δ ξ u i (x)δ ξ u j (x)δ ξ u j (x)dξ. (2.19)
By reinserting the denition of the increments δ ξ u j (x) = u(x + ξ)u(x) and expanding the products, we obtain:

D I (x) = 1 4 ∂ i φ (ξ)u i (x + ξ)u j (x + ξ)u j (x + ξ)dξ -u i (x) ∂ i φ (ξ)u j (x + ξ)u j (x + ξ)dξ + u j (x)u j (x) ∂ i φ (ξ)u i (x + ξ)dξ -2u j (x) ∂ i φ (ξ)u i (x + ξ)u j (x + ξ)dξ + 2u i (x)u j (x) ∂ i φ (ξ)u j (x + ξ)dξ .
(2.20)

The integrals in this last expression can be interpreted as convolutions up to a minus sign, knowing that φ is even so its derivative is odd. The resulting convolutions can then be computed eciently using Fast Fourier Transforms (FFTs).

This makes it particularly convenient in the processing of data from simulations with periodic boundary conditions. The same method can be applied to compute D ν or F using FFTs. The code has been developed in C with OpenMP parallelization and uses the FFTW3 library [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. It has been validated against two other independently developed codes for the computation of the energy transfer term D I .

2.3 Localization of potential singular events using the local energy transfer and turbulent force

These tools have been used in the past to track rare and potentially singular events in turbulent elds. The results discussed in this section are extracted from the thesis of Dennis Kuzzay and Paul Debue. They focus particularly on the In the thesis of Dennis Kuzzay, the measurement were performed using Stereo Particle Image Velocimetry (or SPIV). Using this method allowed him to extract velocity elds with the three components of velocity on a 2D vertical plane passing through the center of the setup. The most interesting results for us were obtained from an experiment using a mixture of water and glycerol. This mixture allows to reduce the Reynolds number down to Re = 2πR 2 f /ν ≈ 6000, which allows measurements down to the Kolmogorov scale. Using these measurements, 33 potentially singular events have been found using a criterion based on local energy transfers. Among those, four dierent categories were observed. Visualizations are provided in gure 2.2. Those categories, based on an interpretation using only a 2D slice, are "front","vortex", "jet" and "cusp". The "front" singularities, which include events where two volumes of uid with very distinct velocity appear to interact, are the most common. The evolution of the local energy transfer as a function of scale provided in inset illustrates the behavior discussed previously:

the local energy transfer is still increasing at a scale close to Kolmogorov.

The same events can also be detected using a method based on the turbulent force instead. Among the events detected with extreme value of the turbulent force in norm, we recognize the same 4 categories found by Dennis Kuzzay, as illustrated in gure 2.3. The evolution of the turbulent force as a function of scale is represented in the insets in the same gure. In all the cases illustrated, we observe a rise of the turbulent force as the scale decreases, which is then dampened for a scale of order = 5η (with the notable exception of the event in gure 2.3a). This short rise is similar to the blow up at small scale expected for real singularities. In this case, the divergence is interrupted, although it is not clear if it is by viscosity or if the limited resolution introduces errors in the computation of the turbulent force. This is what is typically expected for quasi singularities, i.e. a singular behavior observed close to the Kolmogorov scale but interrupted either by viscosity or by the resolution of the dataset.

The comparison of the events found using the turbulent force and the local energy transfers shows that we encounter the same types of event. Moreover, the computational cost associated to turbulent forces is higher than for local energy transfer. For these reasons, we prefer using only the local energy transfer in the following.

This study and classication has been expanded to events found in 3D measurements in Paul Debue's thesis. Those results will be discussed in more details in Chapter 4.

The Velocity Gradient Tensor Invariants

Another tool used by Paul Debue in order to characterize the singular events is the Velocity Gradient Tensor (VGT) invariant method. This method is rst described in a paper by Chong [13] and is used to classify the dierent possible topology for the streamlines of velocity. This classication relies on the eigenvalues of the VGT:

(∇u) ij = ∂ j u i . Depending on the degeneracy of these eigenvalues, whether they are complex or real, positive or negative, there are only a nite amount of possible congurations. Assuming those eigenvalues are λ 1 ,λ 2 and λ 3 in no particular order, the VGT invariants P , Q and R are computed as follow:

27R 2 + 4Q 2 = 0 R Q Vortex Compressing Vortex Stretching Sheet Filament
P = λ 1 + λ 2 + λ 3 , Q = λ 1 λ 2 + λ 1 λ 3 + λ 2 λ 3 , R = -λ 1 λ 2 λ 3 (2.21) (2.22) (2.23)
such that the characteristic polynomial is λ 3 + P λ 2 + Qλ + R.

Furthermore, the incompressibility condition translates to P = 0. As a result, the computation of the invariants Q and R simplies to:

Q = - 1 2 Tr(∇u) 2 , R = -det(∇u). (2.24)
(2.25)

With these constraints, only 4 main topologies remain possible. Those are placed on the QR plane in gure 2.4. In detail, they are the following:

The lament: As 27R 2 +4Q 2 < 0, all three eigenvalues are real. With R < 0, only one of them is positive. The uid is compressed in two directions and stretched in a third one, without rotation.

The sheet: Like for the laments, all eigenvalues are real, but as R > 0, two of them are positive. The uid is compressed in one direction and stretched in another, still without rotation.

24CHAPTER 2. Tracking singularities using local energy transfer and the turbulent force

The vortex stretching: As 27R 2 + 4Q 2 > 0, there are two complex conjugate eigenvalues. As R < 0, the remaining eigenvalue is positive. As the name suggest, the corresponding structure is a vortex spiraling inward stretched in the third direction.

The vortex compressing: Like for the vortex stretching, two eigenvalues are complex conjugate. As R > 0, the third eigenvalue is negative. This corresponds to a vortex spiraling outward and compressed in the third direction.

In the same gure 2.4 are drawn three Vieillefosse lines, dened by:

27R 2 + 4Q 3 = cst.
(2.26)

According to a paper by Vieillefosse [START_REF] Vieillefosse | Internal motion of a small element of uid in an inviscid ow[END_REF], these lines described the possible trajectories of an elementary uid particle in an inviscid liquid not subject to anisotropic pressure eects. The corresponding equations would be:

Q = -3R/2, Ṙ = Q 2 /3. (2.27)
(2.28)

In consequence, should the elementary particle of uid follow these trajectory, this would lead to Q → -∞ and R → +∞ in nite time. In other words, this would lead to a singularity. As such, it is pertinent to verify the behavior of the criteria for the detection of singularity on the QR plane. This work has been done in Paul Debue's thesis, by taking conditional average of D I for given Q and R.

The results feature a higher average local energy transfer below the right branch of the Vieillefosse line drawn in red in gure 2.4, in the "sheet" region. Those results are reproduced in this thesis in Chapter 4.

Chapter 3

Measure of local regularity

The quasi singularities can be identied using the local energy transfers as a criterion, as described in the previous chapter. However, such method only provides an incomplete picture of the intensity of the singularities. In particular, while the criterion based on local energy transfers allows to detect intense events, it does not provide a quantitative measurement of how singular those events are. The ideal tool for this task would be a local scalar that gives a measure of the strength of singularities.

This chapter presents the method developed during this thesis to compute such a scalar. In order to reach this method, we rst describe the multifractal framework which provides statistical tools to quantify the global regularity of a function. We then present a more recent concept of active volumes, which is a mathematical tool developed to locate regions of interest in a turbulent velocity eld. These two preexisting concepts are then combined in order to obtain local information about the regularity of velocity elds.

Overview of the multifractal framework

The Hölder exponents or a measure of regularity

In order to quantify the regularity of the velocity eld when searching for quasi singularities, we need a mathematical tool that quanties the regularity of a function. The concept of Hölder exponents (or Hölder coecients) is such a tool. This scalar is a local measure of the regularity of a function. The Hölder coecient of the velocity eld u at point x 0 is dened [START_REF] Mallat | Characterization of signals from multiscale edges[END_REF][START_REF] Mallat | Singularity detection and processing with wavelets[END_REF]4] as the larger value of h(x 0 ) such that there exists a polynomial P n of degree n < h(x 0 ) and a constant C > 0 25 such that for x in the vicinity of x 0 :

|u(x) -u(x 0 ) -P n (x -x 0 )| C |x -x 0 | h(x 0 ) . (3.1)
If the velocity eld is n times dierentiable in x 0 , its Taylor expansion around x 0 can be used for P n which leads to h(x 0 ) > n. Likewise, for 0 < h < 1, the equation 3.1 becomes:

|u(x) -u(x 0 )| C |x -x 0 | h(x 0 ) . (3.2)
When interpreting the equation (3.2), one can observe that it is equivalent to the simple continuity in the limit h → 0 + . Likewise, the limit h → 1 gives the property of Lipschitzianity, which is close to dierentiability. As a result, the Hölder exponent for 0 < h < 1 is a scalar that makes the bridge between these two concepts and gives a measure of the regularity of the velocity eld over a continuous spectrum. This observation can be generalized for h > 1 using the equation (3.1).

Hence, the general interpretation is that the higher the Hölder exponent is, the more regular is the velocity eld. Computing a local Hölder exponent is required to both locate quasi singularities and estimate their strength.

One pitfall however is that this measure of regularity is only able to detect algebraic singularities (a.k.a. cusps) but not oscillating singularities (a.k.a. chirps).

Those two types of singularities are illustrated in gure 3.1. For oscillating singularities, the Hölder exponent only evaluates the regularity of the envelop. As a result, the measure of the regularity of a chirp event using the Hölder exponent gives an inaccurate estimation of its strength, in particular with regard to the strong associated gradients. A generalized measure including oscillating singularities exists (see e.g. [START_REF] Jaard | Wavelet methods for pointwise regularity and local oscillations of functions[END_REF]1]) but will not be considered during this study. This will eectively limit the reach of all our tools based on Hölder regularity to cusps.

Ideally, we hope that if the INSE generate chirp singularities, their corresponding quasi singularities observable in numerical and experimental datasets would still be detectable using the tool for cusps.

The multifractal framework

Hölder exponents were rst used in turbulence in the context of the multifractal framework. This model has been introduced by Parisi and Frisch [START_REF] Parisi | On the singularity structure of fully developed turbulence[END_REF] and is based on the power law behaviour of the structure functions [START_REF] Frisch | Turbulence: the legacy of AN Kolmogorov[END_REF][START_REF] Monin | Statistical uid mechanics[END_REF]:

δu p ∝ ζp . (3.3)
The Hölder exponent from equation (3.2) is interpreted as a measure of the local scale invariance of the velocity increment δu (x) ∝ h(x) . The relative distribution of the Hölder exponents in space is represented by the multifractal spectrum or singularity spectrum. Given the subset of Hölder exponent ∆(h) dened by ∆(h) = {x| δu (x) ∝ h }, (3.4) the multifractal spectrum D(h) is the Hausdor dimension of ∆(h) :

D(h) = Dim H (∆(h)). (3.5)
The Hausdor dimension is one avor of fractal dimension. Where a line has a dimension of 1 and a surface a dimension of 2, fractal dimensions are a generalization of this concept for more complex structures that do not fall in any simple category. A more easily interpreted concept, used further in this chapter, is the box counting dimension.

Another interpretation can be made in terms of "frequency". In the case of the subsets denes in equation (3.4), the closer the value D(h) is to the dimension of the space, the more ubiquitous are the points with Hölder exponent h.

From here, the rate function of the Hölder exponent is derived as:

Prob [ln( δu ) = h ln( /L)] ∼ e ln( /L)(D-D(h)) = L D-D(h) . (3.6)
where D is the dimension of the full space. By injecting equation (3.6) in the computation of the average in equation (3.3), we can demonstrate that the multifractal spectrum and the scaling exponents of the velocity structure function are connected through a Legendre transform [START_REF] Frisch | Turbulence: the legacy of AN Kolmogorov[END_REF]:

ζ p = min h (ph + D -D(h)) (3.7)
As pointed out by Muzy et al. [START_REF] Muzy | Multifractal formalism for fractal signals: The structure-function approach versus the wavelet-transform modulusmaxima method[END_REF], multiple diculties arise in derivating the multifractal spectrum directly from the structure function exponents. In particular, the potential existence of regular regions with h > 1 would mean using the denition from equation (3.1) instead of equation (3.2). This would impact the computation of the coecients ζ p and thus the derivation of the multifractal spectrum through a Legendre transform [6,[START_REF] Muzy | The multifractal formalism revisited with wavelets[END_REF]3].

Wavelet transforms

The velocity increments used until now can also be interpreted as a particular type of wavelet transforms [START_REF] Grossmann | Decomposition of functions into wavelets of constant shape, and related transforms[END_REF][START_REF] Meyer | Ondelettes[END_REF][START_REF] Meyer | Wavelets and operators[END_REF]. More explicitly, in 1D, the increment of a signal δ f (x 0 ) = f (x 0 + )f (x 0 ) can be written as:

δ f (x 0 ) = 1 ∆ (1) x -x 0 f (x)dx (3.8)
where ∆ (1) (x) = δ(x + 1)δ(x).

Wavelet transforms were introduced by Morlet and Grossman in the context of seismology as a powerful tool in signal processing. The transformation consists in decomposing a signal using a set of functions characterized by a position x 0 and a scale . This property makes it suitable to study local scale invariance, which is the fundamental property of multifractals.

Hence, Arneodo et al. [3] proposed a new method to apply the multifractal analysis that does no share the problems of the structure function based method and can thus be applied on functions with h > 1. This methods relies on continuous wavelet transforms while the classical method used box increments.

In this thesis, we will only consider two dierent types of wavelets: the Gaussian wavelet and the Mexican wavelet. They are respectively dened using the Gaussian and the mexican hat as smoothing functions:

φ Gauss (x) = e -x 2 /2 , φ M ex (x) = (2 -x 2 )e -x 2 /2 .
(3.9) (3.10) In one dimension, this leads to:

Ψ Gauss (x) = dφ Gauss dx (x) = -xe -x 2 /2 , Ψ M ex (x) = dφ M ex dx (x) = x(x 2 -4)e -x 2 /2 . (3.11) (3.12)
Those two wavelets are represented in gure 3.2.

The wavelet transform is then computed by:

T Ψ {f }(x, ) = 1 Ψ( x -y )f (y)dy. (3.13)
Extending this concept in several dimensions imposes to compute elementary wavelets that are partial derivatives of the smoothing function with respect to each direction:

Ψ i = ∂φ ∂x i . (3.14)
The total wavelet transform is then expressed as a vector. In 3D, this results CHAPTER 3. Measure of local regularity in:

T Ψ {f }(x, ) =    T Ψ 1 {f }(x, ) = -3 Ψ 1 ( -1 (x -y))f (y)dy T Ψ 2 {f }(x, ) = -3 Ψ 2 ( -1 (x -y))f (y)dy T Ψ 3 {f }(x, ) = -3 Ψ 3 ( -1 (x -y))f (y)dy    . (3.15)
As such, multidimensional wavelet transforms capture the variation of a function along all directions. This makes the norm of the wavelet transform an equivalent to the increments.

One particular advantage of the wavelet transforms is the ability to deal with the h > 1 case which was problematic for the structure function method. Indeed, a judicious choice of wavelet can cancel out the polynomial P n from the signal. Formally, let's take a function f with a local Hölder exponent h(x 0 ) ∈]n, n + 1[. Taking the same interpretation as before, this means that around x 0 :

f (x) -f (x 0 ) -P n (x -x 0 ) ∼ |x -x 0 | h(x 0 ) . (3.16)
Now let us take a wavelet Ψ with n Ψ cancellation, i.e. such that:

∀k ∈ 0, n Ψ , x k Ψ(x)dx = 0.

(3.17)

Taking the wavelet transform of f around x 0 then gives [START_REF] Decoster | Analyse multifractale d'images de surfaces rugueuses à l'aide de la transformation en ondelettes[END_REF]4]:

if n Ψ > h(x 0 ) : T Ψ {f }(x) ∼ |x -x 0 | h(x 0 ) , if n Ψ < h(x 0 ) : T Ψ {f }(x) ∼ |x -x 0 | n Ψ . (3.18) (3.19)
Using Gaussian wavelets allows us to get good results for h < 1 and using Mexican wavelets extends this range to h < 3.

The WTMM method

What follows is a brief overview of the Wavelet Transform Modulus Maxima method or WTMM. The reader is invited to peruse the original references for the details and original implementation (see [START_REF] Decoster | Analyse multifractale d'images de surfaces rugueuses à l'aide de la transformation en ondelettes[END_REF] for 1D and 2D WTMM and [START_REF] Kestener | Analyse multifractale 2D et 3D à l'aide de la transformation en ondelettes : application en mammographie et en turbulence développée[END_REF] for the 3D generalization).

The general idea of the 1D method consists in extracting the local Maxima of the Modulus of the Wavelet Transform (hence the name) of the signal f at any given scale . These maxima are linked across scales to form lines of maxima Λ k that converge when → 0 to the position of singularities. The position of the local maximum of the line Λ k at scale is noted Λ k ( ).

The line of maxima corresponding to a singularity of Hölder exponent h < n Ψ veries for small scales:

M Ψ {f }(Λ k ( ), ) ∼ h (3.20)
where M Ψ {f } is the modulus of the wavelet transform T Ψ {f } dened in equation (3.15). The wavelet transform skeleton is then dened as the set of lines of maxima {Λ k }.

The WTMM method is a statistical tool using this skeleton to determine the multifractal spectrum D(h). Calling Λ( ) the set of lines, subset of {Λ k }, that exists at scale , we dene the partition functions:

Z(q, ) = Λ k ∈Λ( ) M Ψ {f }(Λ k ( ), )) q (3.21)
where q ∈ IR.

Those partitions functions behave as a power law for small scales, leading to the denition of the exponents τ (q): Z(q, ) ∝ τ (q) .

(3.22)

The singularity spectrum is then obtained by Legendre transform:

D T (h) = min q (qh -τ (q)). (3.23)
In practice, we bypass the Legendre transform by computing:

h(q, ) = Λ k ∈Λ( ) ln |M Ψ {f }(Λ k ( ), )|W(q, Λ k , ), (3.24)
D(q, ) =

Λ k ∈Λ( ) W(q, Λ k , ) ln W(q, Λ k , ) (3.25)
where

W(q, Λ k , ) = 1 Z(q, ) M Ψ {f }(Λ k ( ), ))
q is a Boltzmann weight computed from the wavelet transform skeleton. These quantities behave with respect to the scale as follows:

h(q, ) ∼ h(q) ln( ),

D(q, ) ∼ D(q) ln( ).

(3.27)

By combining the exponents h(q) and D(q) obtained as functions of the parameter q, we can recover the singularity spectrum D(h). Note that the multiscale nature of turbulence means that this spectrum is not independent of the scale at which the power laws are measured. We do not expect to observe the same multifractal spectrum in the inertial range and in the dissipative range. By simplicity, we omit the dependency in scale unless required.

The generalization of this method for scalar and vector elds in several dimensions follows the same general construction, the main dierence coming from the determination of the skeleton {Λ k }. In particular, applying the method to a 3D velocity eld u requires the use of a tensorial wavelet transform, dened by:

T Ψ {u}(x, ) =    T Ψ 1 {u 1 }(x, ) T Ψ 1 {u 2 }(x, ) T Ψ 1 {u 3 }(x, ) T Ψ 2 {u 1 }(x, ) T Ψ 2 {u 2 }(x, ) T Ψ 2 {u 3 }(x, ) T Ψ 3 {u 1 }(x, ) T Ψ 3 {u 2 }(x, ) T Ψ 3 {u 3 }(x, )    .
(3.28)

The process of determining the skeleton of maxima and the multifractal spectrum requires the Singular Value Decomposition of T Ψ {u}, leading to:

M Ψ {u} = max j σ j (3.29)
where σ j are the singular values of T Ψ {u}. The computation of the multifractal spetrum D(h) can then be performed using this new denition of M Ψ {u} in equations (3.24) and (3.25).

All these methods have been implemented in the program xsmurf [START_REF] Kestener | xsmurf v1.0[END_REF], which we use in the sequel for all computations related to the WTMM method, with only minor modications to the software.

A local estimate of Hölder exponents

The WTMM method is based on statistical properties, which makes it very robust with respect to noise and limited statistics. However, using this method only yields global information about the regularity of the velocity eld, through the multifractal spectrum. Meanwhile, in order to study the extreme irregular structures of the velocity eld, we need local information about the Hölder exponents.

Some methods have been designed to provide such information based on the chains of maxima from the WTMM method (see e.g. [START_REF] Struzik | Local eective Hölder exponent estimation on the wavelet transform maxima tree[END_REF]). The general idea behind those methods is to extract the power law behaviour from a single chain of maxima in order to determine the local Hölder exponent of the corresponding singularity. However, this involves at least partially giving up on the advantage of the statistical method. Furthermore, the method relies on a fractal behavior robust in scale. With its multiscale nature, turbulence does not fulll the conditions for this method. One objective of this thesis being the study of experimental data, we need to preserve the robustness with respect to experimental noise of the statistical method. Moreover, our method must take into account the multiscale nature of turbulence.

Our objective in the following is to devise a local statistical analysis, which would allow us to evaluate the local Hölder continuity while keeping the robustness of the WTMM method [START_REF] Nguyen | Local estimates of Hölder exponents in turbulent vector elds[END_REF]. The idea of the method is to combine the results from the WTMM, which gives information about the distribution of the Hölder exponent in the space, with a local method delimiting regions of interest in space closely related to Hölder exponents. In this process, we eventually come back to an intuitive approach of the local regularity: the larger the local gradients, the least regular the eld. The objective of the method is to transform this intuition into a rigorous, reproducible method, which would associate a value of Hölder exponent in accordance to the intensity of the local gradients or the amplitude of the local velocity increments.

Wavelet velocity increments

In order to devise this new method, on need to modify the type of wavelet based increments used. Indeed, the method originally developed by Kestener and Arneodo [START_REF] Kestener | Generalizing the wavelet-based multifractal formalism to random vector elds: application to three-dimensional turbulence velocity and vorticity data[END_REF] denes velocity increments using the singular values of the matrix |G ij |.

The method that we propose will require to use the local velocity increments.

As a result, the use of the singular values might introduce undesired noise. This motivates us to dene the wavelet velocity increments δW (u) through the smoothed velocity gradient G ij and its symmetric and antisymmetric parts:

G ij = dy 3 ∇ j Ψ y u i (y), S ij = 1 2 G ij + G ji , A ij = 1 2 G ij -G ji , (3.30)
(3.31) (3.32) where the base wavelet Ψ(x) will be chosen later.

The normal, longitudinal and perpendicular (also referred to as "transverse") wavelet velocity increments are dened as:

δ L W (u) = max ij |S ij |, δ P W (u) = max ij |A ij |, δ G W (u) = δ L W (u) 2 + δ P W (u 2 ) 1/2 . (3.33) (3.34) (3.35)
By construction, the longitudinal component δ L W (u) contain the information about the potential part of the velocity eld while the perpendicular component δ P W (u) provides information about the vortical component. Their combination δ G W (u) is thus a good choice to provide general scaling properties of the velocity eld. In the following, the subscript G will be droped when referencing to the global velocity increment δ G W (u).

This denition uses neither derivatives over the velocity eld nor singular values. This makes our wavelet based velocity increment the smoothest choice possible to quantify the scaling properties of both the potential and the rotational part of a velocity eld.

In the WTMM method described in section 3.1.4, the three type of velocity increments dened in equations (3.33)-(3.35) can be used as a substitute of the velocity increment based on singular values. One might wonder whether the wavelet transform skeleton {Λ k } needs to be recomputed using these new increments. By simplicity, we made the choice of keeping the same skeletons as the original method. The results presented in the following of the thesis will show that this choice preserves the multifractal spectra obtained through the standard WTMM method.

Active volumes and multifractal

The modied WTMM only yields global information about the regularity of the velocity eld. In order to reach local information, we need to combine these results with a more local concept. This new source of information comes from the denition of actives volumes for a turbulent velocity eld [START_REF] Nguyen | Local estimates of Hölder exponents in turbulent vector elds[END_REF]. The original motivation for the development of active volumes is reported in annex B of [START_REF] Nguyen | Local estimates of Hölder exponents in turbulent vector elds[END_REF].The method delimits the so called active regions by clever thresholding on the velocity increments. In this aspect, it appears closely related to our objective of estimating the local regularity based on the intuition of strong velocity increments corresponding to more singular regions.

We rst dene a family of thresholds T p,dp as:

T p,dp ( ) = (δW (u)) p+dp (δW (u)) p 1/dp .

(3.36)

The limit as dp → 0 can be computed as:

T p ( ) ≡ T p,0 ( ) = exp ln δW (u) δW (u) p δW (u) p .
(3.37)

On the other hand, assuming the velocity eld u is multifractal, i.e. δW (u) p ∝ ζ(p) , then we have the scaling:

T p,dp ( ) ∝ (ζ(p+dp)-ζ(p))/dp . (3.38)
Hence, the limit dp → 0 gives the scaling:

T p ( ) ∝ h(p) , (3.39)
where h(p) = dζ(p) dp is obtained in a way alternative to the WTMM method. Note that the expression of T p in equation (3.37) is analogous to the expression of equation (3.24) used in the WTMM method, so we do not expect this method to yields dierent results.

The active volume A p is then dened as:

x ∈ A p if f δW (u)(x) c p T p ( ), (3.40)
where c p is a scale independent constant, to be determined later.

Given the scaling law of the thresholds T p from equation (3.39), the property of the points in the active volume A p can be rewritten as δW (u)(x)

C p h(p)
with C p a scale independent constant. This property is analogous to the denition of Hölder exponents from the equation (3.2), with a reversed inequality. This similarity can thus be used as a connection point between the information about the regularity of the eld provided by the WTMM method and the local property of active volumes. As a result, we dene the scalar eld h, coined local singularity exponent and by abuse of language local Hölder exponent, such that all elements x in A p verify h(x) h(p). This means that the frontier set of A p , which is also the isosurface corresponding to the velocity increment δW (u) = c p T p,0 ( ), is the set of points x of local singularity exponent h(x).

The wavelet based velocity increment δW (u) is continuous in space. This property is thus transferred to the singularity exponent h. This property diers from the multifractal framework, and instead would be more accurate to a multifractional framework [5]. This dierence means that the quantity h is not a real measure of the local Hölder exponent. We will however show that it shares some properties of the real Hölder exponent, making it relevant for measuring the local regularity of the velocity eld.

Box counting method

The denition of the scalar eld h(x) is dependent of the choice of the constants c p . In order to properly mesh the local method with the WTMM, we must ensure that the properties of the eld h match the ones of the Hölder exponent as closely as possible. In particular, the "multifractal spectrum" associated with h which dictates its distribution in space must match the results of the WTMM. Expressed using the fractal dimension of the isosurfaces of h, this means:

Dim x| h(x) = h = D(h). (3.41)
As discussed previously, the isosurface of h = h(p) is dened as the frontier of the volume A p . This frontier is an isosurface of the velocity increment δW (u), such that:

Dim ({x|δW (u)(x) = c p T p ( )}) = D(h(p)). (3.42)
In practice, the fractal dimension can be estimated using a box counting method. The principle of this methods is to count the number of boxes of size a required to recover a set S. For boxes small enough, this number follows the scaling N (a) ∝ a -d where d = D bc (S) is the box counting dimension of the set S.

The intuitive interpretation of this method relies on the following. Let us consider a line segment of length b. This is a structure of dimension 1, and it takes b/a boxes of size a in order to cover it, which gives a scaling N (a) ∝ a -1 . Likewise, a square of size b takes b 2 /a 2 boxes of size a to cover. This gives a scaling N (a) ∝ a -2 for a structure of dimension 2.

The same process can be applied to cases of non integer dimension, which are frequent for fractal structures. As such, the box counting algorithm consists in counting the boxes needed to recover a set and express this number as a function of the box size in order to determine the corresponding power law. The gure 3.3 illustrates the process applied to a well known fractal, the Sierpinski triangle [START_REF] Falconer | Fractal geometry: mathematical foundations and applications[END_REF]. Based on just the three steps represented, the measured fractal dimension would be d ≈ 1.6. As the Hausdor dimension of this fractal is actually ln(2)/ ln(3) ≈ 1.58, this is a good example of the accuracy of the method.

In this work, the concrete implementation of the box counting algorithm consist of a xed grid scan: using a grid of step a starting at the origin, the number N (a) is computed by counting the number of boxes covering the measured set [START_REF] Mandelbrot | The fractal geometry of nature[END_REF]. The mathematical denition of the box counting dimension would normally require us to determine an optimal covering, i.e. covering the structure using the minimal number of boxes, including rotations and osets. Using a xed grid might impact the precision of the measured dimension. However, this choice has the advantage of a very fast implementation, so we rely on the large amount of statistics in order to compensate for the inaccurate count of the optimal number of boxes.

Hence, applying the box counting to compute the dimension of the left hand side of equation 3.42 leads to:

D bc (c p T p ( )) = D(h(p)).

Statistical tting of the constants

Most of the terms in the equation 3.43 can be computed. Indeed, the term T p ( ) can be computed from the velocity eld using equation (3.37), the multifractal spectrum D(h) is accessible through the WTMM method and the function D bc (T ) can be evaluated by using box counting. The term h(p) is accessible both from the WTMM and the power laws on the T p ( ) as described in equation (3.39). We can choose to use the scaling on the T p ( ) for simplicity, as the two methods should be equivalent.

As a result, we can use equation (3.43) to solve for c p . As the functions involved are not monotonous, there can be several solutions. Noting T max the isovalue of velocity increments corresponding to the largest box counting dimension, we lift the ambiguity by imposing c p T p ( ) < T max for p > 0 and T max and c p T p ( ) > T max for p < 0. This choice originates from the fact that a larger velocity increment should correspond to a less regular point of the velocity eld. Matching the dimensions from equation (3.43) gives a relationship in the form of:

h (p) = f (c p T p ( )). (3.44)
The dependency with respect to of the Hölder exponent h(p) and of the function resulting of the matching of the fractal dimensions are made explicit using a subscript. This provides an unique denition of c p .

It is worth noting that the computation of the coecients c p is not technically required to recover the scalar eld h(x). Indeed, it is possible to derive a relationship between the isovalue T and the local singularity exponent h using the following equation:

D bc (T ) = D( h), (3.45)
which is derived from equation (3.41). As before, we lift the ambiguity by imposing h(T ) < h(p = 0) for T > T max and reciprocally. This leads to the equation:

h (x) = f (δW (u)).
(3.46) using the same matching function f . In summary, the algorithm to evaluate maps of local singularity exponent h at a given scale is performed as follow:

1. First the WTMM method is applied using the velocity increments dened in equation (3.35). This provides the multifractal spectrum of the velocity eld

D(h).

Given that the tting of the power law is performed around the scale , the results are specic to this scale.

2. In parallel, we compute the thresholds T p from equation (3.36) and t their power law behaviour around the same scale using equation (3.39).

3. Finally, we compute the box counting dimensions D bc (T ) for isosurfaces of velocity increment δW (u) = T .

4. The last step is to aggregate the results of those three steps using equation (3.46) to get a function matching a value of h to velocity increments. This directly allows us to compute maps of local singularity exponent from a snapshot of velocity eld.

Note that step 2 is only required in order to compute the coecient c p . While the knowledge of the coecients c p is not required to compute the maps of singularity exponents at a given scale, it is still pertinent. Indeed, these coecients are independent of the scale . This property can be intuited from equation (3.40) as both the velocity increments and the thresholds follow the same power law as a function of the scale. This allows to compute the function f from equation (3.44) at a dierent scale without needing to recompute the multifractal spectrum at this new scale, thus skipping steps 1 and 3. This will be proven to be useful as step 1 is both expensive in computations and sensitive to the tting range used.

The invariance of the coecients c p with respect to scale will be further discussed in the following. Note however that the method exposed here does not provide any guarantee for the behavior of the coecient c p as a function of the Reynolds number.

Scalings in the dissipative range

One of the objective of this thesis is to measure the regularity of the eld in the dissipative range. For this purpose, we would need to apply the tools described in this chapter in a range without the inertial scaling laws. In absence of a robust scaling law, one might doubt about the applicability of the WTMM method described in section 3.1.4. Likewise, the scaling of equation (3.39) is necessary for the local method to work.

However, we believe that the method described here can work even outside of the inertial range. The reason for this comes from an universal scaling law exhibited in a paper by Geneste et al. [START_REF] Geneste | About universality and thermodynamics of turbulence[END_REF]. In terms of structure functions, this universal law for turbulence is expressed as:

β(Re) log S p ( ) S p/3 3 ( ) = H(β(Re) log( /η), p) (3.47)
where the S p are the structure functions, H is some function, and β bears the dependency in Reynolds.

This last function is empirically found to be β = 4/(3 log(R λ )) where R λ is the Taylor based Reynolds number [START_REF] Geneste | About universality and thermodynamics of turbulence[END_REF].

The interesting point about this universal law for our purpose is that it includes a large range over which the extended self similarity holds, i.e. such that: S p ( )

S p/3 3 ( ) ∝ h(p)-ζ(3)p/3 . (3.48)
where the exponent h(p)ζ(3)p/3 is invariant with respect to scale. Such scaling law would be enough for the methods described in this chapter to be functional.

We can rewrite this universal law for our thresholds T p . The universality would be expressed as:

β(Re) log T p ( ) T 0 ( ) = G(β(Re) log( /η), p) (3.49)
where we choose the threshold at p = 0 as a reference and G is some function. This choice is motivated by the fact p = 0 would correspond to the most probable Hölder exponent. With this choice, we would have the same extended tting range as for the structure functions:

T p ( ) T 0 ( ) ∝ h(p)-h(0) . (3.50)
The extended range over which the scaling is veried allows us a better estimation of the power laws, although the nal accuracy of the Hölder exponents is limited by the estimation of h(0). Indeed, when using the method in the dissipative range, this last term has to be computed using a narrow range of scales.

This universal law can also provide some insight about the behavior of the coecients c p with respect to the Reynolds number. Indeed, taking the denition of active volumes in equation (3.40) and dividing by T 0 , we get:

x ∈ A p if f δW (u)(x) T 0 ( ) c p T p ( ) T 0 ( ) . (3.51)
The universal law from equation (3.49) tells us that the right hand side of this inequality behaves as follows:

c p T p ( )/T 0 = c p exp 1 β(Re) G(β(Re) log( /η), p) . (3.52)
At the same time, the left hand side can be interpreted as the ratio between the velocity increment δW (u) and the most probable increment T 0 . Using the same thermodynamic analogy as in [START_REF] Geneste | About universality and thermodynamics of turbulence[END_REF], we would expect this ratio to follow a scaling of the type exp(-∆F/β(Re)) where ∆F is a dierence in free energy.

Hence, for the expression of equation (3.52) to have the correct behavior, we would expect the coecients c p to verify the following for all Reynolds numbers:

c β(Re) p = cst(p). (3.53)
Taking into account that β(Re) ∼ 1/ log(Re) [START_REF] Geneste | About universality and thermodynamics of turbulence[END_REF], this would lead to:

c p = Re α(p) (3.54)
with a function α(p) to be determined.

In this thesis, the datasets analyzed are produced at only two dierent Reynolds numbers. This will be insucient to verify this property with a sucient accuracy, leaving it as a work for the future.

Validation of the method on fractional Brownian motion (fBm)

The method explained above is designed to yield a measure of the local regularity of a velocity eld through the scalar h, assimilated to a local Hölder exponent.

However, while this scalar satises some properties of the true Hölder exponent by construction, we are not able to provide a mathematical proof that it is the Hölder exponent. As a result, a validation is required in order to use the local singularity exponent as a criterion for the detection of singularities in a real velocity eld. In addition to provide a validation of the method, this section is meant as a step by step illustration of the algorithm.

Denition and implementation of fBm

A rst test of our method is to apply it to synthetic elds with well known fractal behaviour. In this case, a simple choice is the case of the fractional Brownian motion (or fBm), which has also been used in [START_REF] Kestener | Generalizing the wavelet-based multifractal formalism to random vector elds: application to three-dimensional turbulence velocity and vorticity data[END_REF]. Such a eld is designed to have an uniform regularity. As such, its fractal behaviour can be described by a uniform Hölder exponent H, called in this case the Hurst exponent [START_REF] Mandelbrot | Fractals: Form, Chance, and Dimension[END_REF].

This means that its multifractal spectrum in dimension d degenerates to a single point:

D(h) =    d if h = H, -∞ otherwise. (3.55)
By Legendre transform, this leads to the following formula for the exponent τ (q):

τ (q) = qH -d. (3.56)
The objective of the following test is to recover the value of H from a generated database of fBm. A similar test was performed to validate the original WTMM method in [START_REF] Kestener | Analyse multifractale 2D et 3D à l'aide de la transformation en ondelettes : application en mammographie et en turbulence développée[END_REF], which justies that we use the same test to validate our variation of it, using the dierent components of the wavelet based velocity increments. The fBm are generated in 2 dimensions with 3 components of velocity (2D3C) with a value of H = 1/3. The generation is performed by ltering the Fourier transform of a white noise, which gives a good approximation of the expected scale invariance properties.

Test of the global MFR analysis

A total of one hundred 2D3C fBm are generated in squares of size unity with a resolution of 4096 × 4096. The output of our method using the velocity increments dened in equations (3.33) to (3.30), as well as the original WTMM method using Singular Value Decomposition are shown in gure 3.4.

The results from our method match the theoretical prediction, although with less accuracy than the original method. Indeed, the value of Hölder exponent found using the new method is H ≈ 0.31, while the classic WTMM gives a result of H ≈ 0.33. While it is expected that the multifractal spectrum computed from a nite number of snapshot does not fully degenerate as in equation (3.56), we also observe that the spectrum computed using our method spreads over a larger range of Hölder exponents. The loss of accuracy might be partially explained by the lack of the third dimension hindering the computation of the perpendicular component from equation (3.34).

Test of the local estimate

The next step consists in extracting the local singularity exponent h from the fBm snapshots. For these computations, we set a scale = 2.8 × 10 -3 corresponding to 11.4 grid points for the resolution of the eld. This guarantees that the wavelets are well resolved and allows to probe small structures. The dierent quantities to compute in order to apply the method are the box counting dimensions associated with the isosurfaces of velocity increments D bc (T ) (equation (3.45)) shown in gure 3.5 and the thresholds T p ( ) (equation (3.37)) shown in gure 3.6. Both are computed using the full statistics. We check that the power law behavior of the threshold from equation (3.39) is consistent over a large range of scales. Furthermore, we observe little dependency of h with respect to p as expected for a monofractal, with a measured value of h ≈ 0.34. This is close to the expected value of H = 1/3.

At this point, the method described in section 3.2.3 would involve matching the multifractal spectrum from gure 3.4 with the box counting dimension from gure 3.5. However, as the spectrum is degenerated in this case, this step is impractical. A substitute for this monofractal case would be to choose the c p such that for all p, c p T p ( ) = T max , eectively guaranteeing the degeneracy of the spectrum. Furthermore, we take the freedom to renormalize the box counting dimension D bc (T ) such that D bc (T max ) = 2. This is justied by the fact that getting measurements of box counting dimension that exactly reach 2, the dimension of the space, is very dicult given the accuracy of the power laws tted. If the fractal dimensions from the WTMM and the box counting method are not properly renormalized, this would lead to an unnatural jump in the matching function f obtained through equation (3.46).

The monofractal nature of fBm makes maps of h of little interest. A more instructive test is to verify the scale invariance of the coecients c p in order to compute elds of local exponent h at a dierent scale. For this purpose, we compute the box counting dimension at the lower scale = 8.2×10 -4 corresponding to approximately 3.4 grid points. The objective is to verify that the equation (3.43) The box counting dimensions of the isosurfaces of "velocity" increments are shown in gure 3.7. By combining it with the c p computed at the scale and the scaling exponents from the T p ( ), we can recover the multifractal spectrum using equation (3.43). The results are displayed in gure 3.8. They match the theoretical spectrum, with a slightly higher value of Hölder exponent than the value found by the WTMM method. This is coherent with the previous results and provides a rst validation of the scale invariance of the coecients c p .

10 -3 10 -2 10 -1 ℓ 10 -1 10 0 T p (ℓ) p p = -0.6 p = -0.4 p = -0.2 p = 0.0 p = 0.2 p = 0.4 p = 0.6 p = 0.8

Summary of this chapter

In this chapter, we introduced the multifractal framework and the WTMM method used to compute multifractal spectra [START_REF] Kestener | Analyse multifractale 2D et 3D à l'aide de la transformation en ondelettes : application en mammographie et en turbulence développée[END_REF]. These can provide global statistical information about the regularity of a velocity eld. We combined these concepts with the notion of active volumes, providing local information about the turbulent velocity elds. The resulting method allows us to estimate the local Hölder exponents of a velocity eld, provided statistical information about the dataset.

The validity of this method has then been tested using a synthetic case.

Chapter 4

Applications to turbulent elds

In this part, we apply the methods described previously on turbulent velocity elds originating from Direct Numerical Simulations (DNS).

The rst objective is to further validate the method based on local Hölder exponent. For this aim, it is ideal to use elds for which we can precisely control the parameters, hence the numerical simulation. The main advantages to DNS elds in this case are the absence of noise and the well controlled and periodic boundary conditions. For the sake of the validation of the method, we can use the turbulent eld in the inertial range, as this provides a range of scale over which we know that the Hölder regularity is of the order of h = 1/3 (see [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous uid for very large Reynolds numbers[END_REF]).

Once the validation is done, the second objective is to detect and characterize potential singularities using both the local Hölder regularity and the local energy transfers. For this purpose, we need to probe velocity elds in the dissipative range, down to the Kolmogorov scale. This study will likewise be performed on numerical velocity elds, but can be extended to experimental results.

Simulation in the inertial range 4.1.1 Description of the datasets

The turbulent velocity elds are generated using two DNS in a periodic box of size (2π) 3 . The parameters for the simulations are presented in Table 4.1. The code used is a pseudo spectral code adapted [9] with a Taylor-Green forcing. The code is parallelized using OpenMP and the Fourier Transform are computing using the FFTW library [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. Such forcing is chosen in order to make the parallel between the simulations and the experimental elds that have been extracted from a Von Kármán experiment. The experimental results which have been mentioned in Chapter 2 will be the object of further comparisons in this chapter.

The averaged energy spectra corresponding to those two simulations are presented in gure 4.2. This illustrates that the rst simulation (Run I) features a large range of scale in the inertial regime while the second simulation (Run II) is resolved beyond the Kolmogorov scale. As such, we can use the rst simulation in order to validate the method using a well known power law behavior. Meanwhile, the extensive dissipative range allows us to probe the actual singular behavior of the ow.

The simulations have been run on the Jean-Zay supercomputer, equipped with Intel Cascade Lake 6248 processors. Using 20 cores for each simulation, the total integration time is of the order of 20,000 CPUh.

In this section, we focus on the results from the Run I. The simulation has been integrated over a total time of 40s, or more than 400 eddy turnover times t η = (ν/ ) 1/2 . As can be observed from the plot of the time evolution of energy dissipation represented gure 4.1, this integration time includes a transitory period of about 10s. The statistics used for the inertial case are computed from 10 snapshots of velocity elds extracted evenly from the end of the transitory period.

As such, two successive elds are separated by a period of about 3s, or more than

Validation of the WTMM method in the inertial range

The rst step for the validation of the method in the inertial range consists in applying the modied WTMM method to the velocity elds from Run I. This is both a verication that the modication we made to the method does not compromise its accuracy in the multifractal case, and a step towards the application of the local method. Furthermore, we can take the opportunity to test the multifractal behavior of the dierent components of the velocity increments (i.e the longitudinal and perpendicular components) independently.

The method is applied on 10 snapshots of velocity eld decorrelated in time.

Given the resolution of 768 3 , this represent more than 4 billions data points, enough to guarantee the convergence of results presented in gure 4. The multifractal spectrum computed is centered around a value of h = 0.40 for all components but the perpendicular one. This is higher than the value of h = 1/3 expected from the Kolmogorov 41 theory [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous uid for very large Reynolds numbers[END_REF]. We however observe a very close agreement between our modied method using the global increment δ G W (u) and the unmodied version of the WTMM, which is enough to validate our modication. The discrepancy of the results compared to the expected singularity exponent can be attributed to the anisotropy of the ow chosen. Indeed, using the same method, Arneodo [2] has found a most probable Hölder exponent of around h = 0.34 for isotropic turbulent velocity elds. The dierent behavior of the multifractal spectrum for the perpendicular component at this scale may also be attributed to the anisotropy caused by the TaylorGreen forcing.

Computation of the local singularity exponent

Another point of our method that requires validation is the local estimate of the regularity. In particular, we have yet to justify how the scalar eld h dened in the previous part can be used as a pertinent measure of local regularity even if it is not the real Hölder exponent of the multifractal framework. For the purpose of this verication, we rst compute it in the inertial range.

The computation follows the same steps as shown for the fractional Brownian motion in part 3.3. In this case, the scale is set at the value ≈ 0.24, well within the inertial range. The method involves computing box counting dimensions over isosurfaces of the velocity increments. As described in the previous chapter, this process consists in the counting of "boxes" of given size required to recover an entire isosurface. The process is repeated over all snapshots of velocity in the dataset. The total box counts in this case are reported in gure 4.4a for six dierent isovalues. We can verify that these counts follow convincing power laws as a function of the size of the boxes. The exponent of the power laws is the fractal dimension associated to the given isosurface.

The results of the box counting, shown gure 4.4b, can be interpreted in terms of distribution of the velocity increments within space. As such, extreme values of velocity increments, either low or high, are concentrated in smaller structures, leading to a fractal dimension lower than 3. Compared to the box counting results for a fBm from gure 3.5, we observe a range of velocity increments for which the corresponding box counting dimension exhibit a plateau close to 3. As observed in the box counts in gure 4.4a, the corresponding isosurfaces display very similar behavior, with close box counts. This can be interpreted as a approximate way to observe intermittency, with large change of the velocity increments over very short period. Such change would lead to the isosurfaces of increments being close to each other hence the similar box counts. It is however unclear how pertinent this metric alone can be for the study of intermittency. The thresholds T p ( ) are computed over the same dataset. We verify in gure 4.5a that once renormalized by T 0 ( ), the results give a convincing power law over a wide range of scales, compatible with the prediction of equation (3.50). In order to compute the power law exponents, the t are performed rst over the renormalized thresholds, which yields a quantity γ(p) = h(p)h(0). The reference value h(0) is then evaluated by tting a power law over the thresholds T 0 ( ) in gure 4.5b over the inertial range. As a reminder, the thresholds and the corresponding power laws are only used to compute the coecients c p , which are not strictly required to recover the local Hölder exponents. We will however use these results to verify the invariance of c p with respect to scale.

As the multifractal spectrum in this case is non degenerated, we can proceed dierently from the fBm case. We use equation (3.45) to match the fractal dimensions from the WTMM method (gure 4.3) with the ones obtained by box counting (gure 4.4b). This gives us a direct relationship between the velocity increments and the local Hölder exponent h. In the process, the lower values of h must be associated with the larger values of velocity increments, as it represents a more singular velocity eld. The resulting matching function f is represented gure 4.6. The perpendicular and longitudinal increments are no longer treated separately as we are interested in the regularity properties of the full velocity eld.

This process allows us to produce maps of h from the corresponding maps of velocity increments, with one pitfall: the range of velocity increment accessible to this method is limited. In particular, for the boxcounting method to work for an isosurface at value T , a minimum amount of data points with such value of velocity increment is required. Likewise, the WTMM method is also constrained by the amount of statistics available. A lower fractal dimension corresponds to rarer events, which makes the collect of sucient statistics harder. As a result, no reliable fractal dimension has been computed for isosurfaces corresponding to a value T > 0.5, which prevents us from measuring values of h < 0.1. In consequence, the maps of singularity exponent that we can compute have regions where we only know that h < 0.1 without being able to give any precise value.

Comparison with the local energy transfers

The production of these maps of h is the opportunity to compare this measure of regularity to the local energy transfers. As mentioned in Chapter 2, the Duchon Robert term D I of the energy budget is expected to be higher in the regions where the velocity eld is singular. Conrming that this term is correlated to the eld h would validate the use of the latter as a criterion to detect pseudo singularities.

Computing maps of those two quantities in the snapshots from Run I shows that these two quantities are visually correlated, with structures of low regularity measured by h being observed in the same location as structures of extreme energy transfer D I . This is illustrated in gure 4.7. In particular, regions with low h appear to be surrounded by a pair of regions of high energy transfer. This also means that the extrema of those two quantities do not occur exactly in the same location.

Nevertheless, it is still relevant to compute the joint Probability Density Function (PDF) of these two terms in order to conrm their correlation. Such PDF, computed over 50 velocity elds extracted at regular interval over approximately In order to do this, one must account for the errors generated by tting power laws. Indeed, the computation of the multifractal spectrum using the WTMM requires to determine power laws according to equations (3.26) 4.9: Pseudo spectra reconstituted using the coecient c p computed at scale ≈ 0.24 and the boxcounting dimensions of the velocity increments computed at scales ≈ 0.12 (green curve), ≈ 0.24 (red curve) and ≈ 0.46 (blue curve). The scale are also materialized by vertical lines of the same colour on the gure 4.5a.

The coecients c p do not have an obvious physical meaning, making it dicult to determine what range of error is acceptable. As a result, we use a similar method as for the fBm case to test their scale independence. This involves computing the c p at a scale ≈ 0.24 and using the results to infer multifractal spectra of the velocity eld at a dierent scale. Using the coecients c p and the box-counting dimensions of the isosurfaces of velocity increments computed at two dierent scales, we use the equation (3.42) to recover pseudo multifractal spectra at two dierent scales.

The results are presented in gure 4.9 using the scales materialized in gure 4.5a.

The blue curve shows that using a larger scale in the inertial range results in a similar multifractal spectrum. The shift towards Hölder exponents even lower than the value h = 1/3 can be attributed to the proximity to the injection scale.

Conversely, the green curve demonstrate that using a scale at the lower bound of the inertial range shifts the spectrum, with much higher values of h. This is a result of the velocity eld starting to regularize through the eects of viscosity.

Although a direct comparison with the output of the WTMM method would be ideal, too many statistics would be required in order to have accurate results with very small tting range. Considering the qualitative aspect of the measure of local regularity through the singularity exponent h, we estimate that this analysis is enough to justify the scale independence of the coecients c p . The result of this analysis is that using coecient c p computed at a dierent scale will provide qualitatively sound results. Since the complete method is dependent on tting scaling exponents, it is reasonable to use the c p for quantitative measurements.

VGT invariants

As we have analyzed the correlations between the Hölder exponents and the local energy transfers, it can be pertinent to reproduce the study by Paul Debue on the VGT invariants [15] as mentioned in Chapter 2. The objective is to verify that we obtain similar results for the location of high energy transfers on the QRplane. Additionally, we can conduct a similar study focused on the local Hölder exponents.

Because most quantities computed in this thesis are ltered using a Fourier transform, we choose to do the same for the VGT invariants. As a result, instead of the standard velocity gradient, we use a mollied gradient dened similarly to equation (3.30):

Q = - 1 2 Tr(A 2 ), R = -det(A ) (4.1) (4.2)
where:

A ,ij = -∇ j Ψ mex (y)u i (x + y)dy. compared to [15], but the qualitative observation are the same. We can verify that we obtain the same common droplet shape for the distribution of the invariants Q and R in gure 4.10a. We also recover similar results for the local energy transfer in gure 4.10b. Indeed, we have the same behavior with strong positive energy transfer along the Vieillefosse line in the lower right quadrant, corresponding to the limit between sheets and vortex compression, as well as in the vortex stretching region. Meanwhile, negative energy transfers are found for high values of Q, both in the vortex stretching and the vortex compressing regions.

The local Hölder exponents however do not follow the same distribution. Overall, we observe lower regularity whenever the invariants Q and R are far from the average, in particular but not only in the vortex region. This is not unexpected, as high values of the VGT invariants require strong velocity gradients, and hence a We also use this part of the dataset to extract singular events in order to identify the typical structures.

The second part of the simulation consists in saving a velocity eld every 30 time steps, or 0.15t η . The purpose of this part is to analyze the time evolution of irregular events. This would be useful in order to identify the mechanism behind the formation of quasi singularities.

Multifractal spectrum in the dissipative range

One idea in order to accomplish this goal is to apply the method again for the new dataset. As such, we start by using the modied WTMM, while choosing the range of tting for the power laws in the dissipative range. The range chosen is materialized in gure 4.2 and spans from around 3η to 9η. This range of scale is large for the purpose of measuring a trend particular to a single scale, as would be suitable in theory for the method. It is however a decent compromise between the localization of the scale and the accuracy of the power law tted. As we do not expect consistent power law behavior of the velocity increments over this large range of scale in the dissipative range, the WTMM method might not work as well as in the dissipative case.

For similar statistics as for the inertial case, we use 10 snapshots from the Run II, decorrelated in time. With this amount of statistics, we are condent as before of the convergence of the results. The multifractal spectra obtained are displayed in gure 4.12. One rst observation in comparison to the inertial case in gure 4.3 is that the spectra computed using the dierent components of the velocity gradient follow approximately the same curve. In other words, we have recovered isotropy at the dissipative scale. This is an expected behavior, as the small scales structures should not be too heavily inuenced by the large scale forcing. This nevertheless is an important property to note in order to generalize the small scale observations of this chapter.

One other striking observation is that the range of Hölder exponents observed at this small scale has much higher values than what was observed at large scale.

Again, this is coherent with the expectations. As mentioned before, the eect of viscosity is to make the velocity eld more regular. At this scale, the most frequent Hölder exponent observed is around h max = 2.7. In eect, this means that most of the eld appears to be twice dierentiable at this scale. The spectrum is bounded on the right side by h = 3, which is the maximum value measurable using a Mexican wavelet, with three cancellations. In other words, if the actual Hölder exponent of the eld was any higher, we would still only measure a value of h = 3.

On the other hand, the fractal dimension associated with a Hölder exponent of h = 2 is non zero. This means that there are locations in the velocity elds with a lower regularity at this scale.

One other observation about this spectrum is that the fractal dimension does not reach the dimension of the space. We can interpret it as the lines of maxima dened in the section 3.1.4 are not space lling at this scale. This means that the singular events followed across scales by the WTMM method are more sparse. This aspect is detrimental to our method as we need to match the dimension obtained from the box counting algorithm with the present spectra.

Calibration of the local Hölder for the dissipative range

The behavior of the multifractal spectrum at the dissipative scale makes it impossible to apply the method directly. In particular, the fractal dimension computed is much lower than 3, while we expect the box counting algorithm to give fractal dimensions close to 3 for the most common values of velocity increments. As a result, it is not suitable to use the full method described in part 3.2.

In order to overcome this diculty, we use the property that the coecients c p are independent of the scale. The idea is thus to apply our method in the inertial range, where we know the method works properly. From the results we compute the coecients c p . Those coecients can then be used to recover the matching function between the velocity increments at the dissipative scale and the Hölder exponents at the same scale using equation (3.44).

We thus compute the multifractal spectrum and the boxcounting dimensions at the scale ≈ 1.4λ ≈ 19η. The results of these two operations, as well as the Inset: Exponents τ as a function of q.

resulting matching function between the velocity increments and the local Hölder exponents, are shown gure 4.13. Because the inertial range is much shorter than for the previous simulation and the scale at which the computation is performed is much closer to the injection scale, we encounter diculties for this step. In particular, it is dicult to get the WTMM method to converge, as there are much fewer lines of maxima compared to the previous case. In order to help with the convergence, we use all the 40 velocity elds uncorrelated in time from the rst part of the simulation. We further verify that adding 60 more time resolved velocity elds from the second part of the simulation do not modify the results signicantly. This does not prevent the eect of the proximity to the injection scale, which is made apparent by the fact the most probable Hölder coecient in gure 4.13a drops to h = 0.3, compared to the results in the previous simulation. As such, we need to keep in mind for the future interpretation that there may be inaccuracies in the estimation of the local Hölder exponents.

From this matching function, the next step is to determine the coecents c p . In order to do this, we compute the thresholds T p ( ) for this dataset. Once again, the thresholds are renormalized by the reference value T 0 ( ). The results are show in for the dissipative scale. While the inertial range is narrow, the power law over the dissipative range is sustained over about a decade. This justies the use of multifractal based tools for the dissipative range.

Using the same method as for the inertial case, we can derive the values of h(p) from the power laws tted in 4.14 for both = 1.4λ in the inertial range and = 1.8η in the dissipative range. Using equation (3.44) for the rst scale allows us to compute the coecients c p from the thresholds and the matching function of gure 4.13. We can then use the same equation at = 1.8η to recover the matching function in the dissipative range, shown in gure 4.15. The corresponding Hölder exponents are as expected much higher than those computed in the inertial range from gure 4.13c. This is the eect of viscosity regularizing the velocity elds at this scale. In terms of rarity of an event, one must understand a local Hölder exponent of 1.4 at this scale to be as exceptional as a Hölder exponent of 0.2 in the inertial range.

Note that the corresponding local Hölder exponents are lower than the one determined using WTMM in gure 4.12. We do not have a clear explanation for this discrepancy, as for both cases a Mexican wavelet is used, which is supposed to only cause problems when the regularity reaches h = 3. For the purpose of our analysis of the local regularity, we will trust the results from gure 4.15, at least for qualitative interpretations. Even taking into account the uncertainty on the exact values of Hölder exponent, it is clear that there are events with h < 2 at this scale, i.e. such that the velocity eld does not appear to be twice dierentiable.

In the following, we will try to probe for very irregular events at the dissipative scale. In order to give a value of Hölder exponent even if the velocity increment is greater than 0.2, we choose to extrapolate the matching function as shown in gure 4. 15. Events with such low regularity are expected to be rare, and we expect the local Hölder exponent to reach such low value only in very localized regions, as it was the case for the maps shown in gure 4.7. As such, we deem that a simple linear extrapolation based on the last two datapoints should be enough for our purpose. In order to associate a local Hölder exponent to these quantities without the use of extrapolation, one would need to compute thresholds and multifractal spectrum for higher orders of p. However, the amount of statistics available to us is not enough to get converged results for p > 5.

Additionally, we are also unable to map a Hölder exponent to the very small values of velocity increments. This is because such values are also very rare, corresponding to the p < 0 branch of the multifractal spectrum. While this has no eect on our study of very singular events, one can expect the Hölder exponent to

Comparison to the local energy transfers in the dissipative range

Now that we have access to the local Hölder exponents at the dissipative scale, it is important to perform another comparison with the local energy transfer.

Indeed, the purpose of the study performed in section 4.1.4 was to validate the use of the local Hölder exponents computed using this method as a criterion for the detection of quasi singularities. While a good correlation between these two scalars computed in the inertial range is a rst step for validation, the detection of irregular event is to be performed at the dissipative scale.

As a result, we compute the joint PDF of these two terms computed at = 1.8η. location as the maxima of D I . This contradiction is lifted by plotting the same PDF restricted to very low values of Hölder exponents. The results are shown in gure 4.17. We observe that for very small values of Hölder exponent, the correlation is reversed. This shows that, particularly at the dissipative scale, the structures of very low Hölder exponents and very large local energy transfer are not superposed. This observation will be further discussed in the following.

VGT invariant in the dissipative range

When analyzing extreme events at the dissipative scale, we will make use of the VGT invariants in order to get an idea of the topologies most representative of quasi singularities. For this purpose, we need to know the global, unconditioned PDF of the invariants Q and R as a reference. As for the Run I, we use VGT invariants computed over the mollied velocity gradient, as dened in equation From the results, we select all events where the minimum value of h gets below 1.3. This value is chosen to strike a balance between the strength of the events and the number of events found. With such a low value of h at the dissipative scale, only inferred through the matching function of gure 4.15 by using extrapolation, we target very irregular events. On the other hand, this threshold is high enough to allow us to extract 213 distinct events fullling this condition. This will be enough to perform statistics.

We can verify that the events detected are spread evenly among the snapshots of velocity elds, as shown in the histogram of gure 4.19. Likewise, we also verify that the events are not concentrated in a single region of the simulation domain.

This guarantees that there is no bias caused by the forcing chosen.

Before analyzing the events in detail, we can use the VGT invariants to get a rst grasp of the topology of the singular events. We dene the location of the minimum of local Hölder exponent as the center of the event. The VGT invariants Q and R computed at the center of each using the mollied velocity gradient (see equation (4.2)) are represented in gure 4.20. The rst observation is that, with only two exceptions, all events appear to be either in the "vortex compressing" or the "vortex stretching" regions of the QR-plane. Moreover, we observe a bias toward the later as more than 75% of the events detected verify R < 0. We conrm that this bias is stronger than what is observed in the whole domain in gure 4.18. This trend is in agreement with several studies relating singularities to vortex stretching (e.g. [START_REF] Kerr | Evidence for a singularity of the three-dimensional, incompressible euler equations[END_REF] in the case of Euler equations).

Visualization of singular events

The next step in the analysis is to observe the events directly. In order to reduce the inuence of the large scale motion in the visualizations, we subtract the mean velocity over a cubic sub-domain of side 1.7λ centered on the point of minimal Hölder exponent.

For each events, we compute the vorticity, the inter scale energy transfer D I (from equation (2.8)), the viscous dissipation at this scale D ν (from equation (2.12)) and the mollied velocity gradient tensor invariants Q and R (from equation (4.2)).

A rst overview of the events detected leads to two observations. First of all, we verify the conclusion obtained through the use of the VGT invariants in gure 4.20: the streamlines of velocity for most events are typical of vortices. Moreover, the isosurfaces of local Hölder exponents h are oblong, seemingly oriented in the same direction as the lines of vorticity. The other principal observation is that In order to help with the analysis, it is preferable to have all of them oriented in the same way. The way the events are reoriented is as follow:

As the local Hölder exponent h is the most relevant scalar for our study, we dene the direction of a rst axis by using the rst eigenvector e 1 of the covariance matrix of h. Reorienting the events using this axis, the structures of Hölder exponent must be aligned along the x-axis.

As singular events appear to be vortices, we want all vortices to face the same direction. We thus choose the direction of the rst axis such that ω x > 0 at the center.

We observed that the events are not axisymmetrical from the structures of D I and D ν . This justies the denition of a second axis, orthogonal to the rst one. Using the rst eigenvector e 2 of the covariance matrix of D ν , we reorient the events such that e 1 × e 2 is along the z-axis.

Note that those choices do not impose a direction to the second axis, such that rotating the event by an angle π around e 1 still fullls the conditions. Three events reoriented as explained are represented in gures 4.21 to 4.23. Additional events can be found in Appendix A.

Using the data collected, we can pick up the work of classication for the extreme events mentioned in Chapter 2. The thesis by Paul Debue [15] approaches the same problem using extreme events of D I . As was discussed before, those are closely related to the extreme events of Hölder exponents studied here. In that work, two types of vortices are distinguished based on the behavior of the helicity:

the roll vortices and the screw vortices. The roll vortices corresponds to events where the helicity changes sign over the structures. On the other hand, the screw vortices are vortices without change of helicity in the structure. It is mentioned in this work that the dierence between the two types of events might just be a simple Galilean transformation.

In our dataset, the predominant structure is the one corresponding to roll vortices. Such a structure is visible in the streamlines of gure 4.22. In some cases, the change of helicity occurs o centered, away from the location of minimum regularity. This can be attributed to a remnant contribution of the large scale velocity.

Meanwhile, screw vortices rarely appear. One could interpret the aforementioned o centered roll vortices, like for example 4.23, as screw vortices. However, we observe no qualitative dierence between those two types of vortices in any of the observable used here. The conclusion is that, as [15] supposed, those two structures are the same up to a Galilean transformation.

In the work of [15], a third type of structure is identied. The so called "U-turn" is characterized by a sharp change of direction of the velocity streamlines. We do not observe any "U-turns" in our dataset. The closest match are tight vortices, such that they could be seen as "U-turn" if observed at a lower resolution and with the addition of experimental noise. It is however possible that "U-turns" t the criterion based on local energy transfer but not the one based on local Hölder exponents, hence they are not detected in our case. It is also possible that "Uturns" correspond to an earlier phase in the development of a quasi singularity, such that the local Hölder exponent is too high to be detected by our method. Finally, it could also be a dierent structure appearing at higher Reynolds number. If the last hypothesis is true, we would not be able to study this type of event through the simulations performed in this thesis.

Among the recurrent features of the events observed is an asymmetry of the vortices. This asymmetry is particularly manifest in gure 4. [START_REF] Frigo | The design and implementation of FFTW3[END_REF], where the isosurfaces of vorticity are compressed in a pancake shape. This shape is shared by the isosurfaces of Hölder exponent, and the isosurfaces of D I and D ν appear as pairs of pancakes sandwiching the event. This reminds the behavior observed by Vincent and Meneguzzi [START_REF] Vincent | The dynamics of vorticity tubes in homogeneous turbulence[END_REF]. In this study, the authors remark that regions of high vorticity are pancake like at their apparition, before bending and rolling up. A connection between this process and the energy cascade is also made. This could explain why the event in gure 4.22 features higher values of D I than the others. For this to occur at a scale of = 1.8η might be another hint of a structure carrying energy below the Kolmogorov scale.

Typical event

Having extracted a large amount of extreme irregular events, we can try to average them in order to obtain a "typical event". This is made possible by the fact all events have been reoriented along the same axis. This is also justied by the fact that most events detected share a similar structure of "roll vortex". Exceptions are composed of events happening close to each other, leading to complex streamlines pattern. However, in order to avoid introducing any bias in the average, we do not perform any ltering of the events. As the complex, non categorized events do not have clear common features, we expect them to not contribute too much in the average.

The average is taken over all 213 most singular events extracted. As a reminder, this corresponds to all events with a local Hölder exponent h < 1.3. Visualizations of the averaged event are provided in gure 4.24. Note that the quantities displayed are the averaged values of the exponent h, the energy transfer D I and the energy dissipation D ν . As those scalars are not linear functions of the velocity, they do not correspond to the quantities that would be computed from the averaged velocity eld, materialized by the streamlines.

The pattern of streamlines for the averaged event corresponds to the roll vortices, which emphasizes the prevalent character of this type of structure among extremely singular events. As expected, averaging over many events completely removed the contribution of the large scale velocity, such that the helicity changes sign at the x = 0 plane corresponding to the minimum of Hölder exponent. From this structure one might recognize a Burgers vortex. This type of uniformly stretched vortex has been dened by [10] and it is characterized by its vorticity prole:

ω x (r) = Γ 2πν exp - σr 2 2ν . (4.4)
This prole can be tted over the vorticity of the averaged event in the plane x = 0. This yields Γ = 5τ ν and σ = 0.36. While the t is accurate in the z direction, as shown in gure 4.25a, it fails in the y-direction. Because the x component of the vorticity actually changes sign, this t cannot be xed using the non axisymmetric Burger vortex prole. The gure 4.25b illustrates this asymmetry.

Likewise, the averaged D I and D ν show a very strong non axisymmetry. The invariance by a rotation of π around the x-axis can be explained by the way the events are averaged. However, the presence of two structures of D I and D ν on either side of the main axis is also observed in individual events (see gures 4.21 to 4.23). In order to interpret this asymmetry of the averaged event, we can emit the hypothesis that extreme events result of the interaction between several vortices. In the process of averaging, only the main, most singular vortex was conserved. This hypothesis is supported by the presence of nearby vortices in several of the extracted events as can be seen in the example of gure 4.23. Under this hypothesis, the structures of vorticity could be similar to what is observed in the bridge of a vortex reconnection [START_REF] Yao | On singularity formation via viscous vortex reconnection[END_REF][START_REF] Yao | A physical model of turbulence cascade via vortex reconnection sequence and avalanche[END_REF]. In order to rene this hypothesis, we can use time resolved data. Indeed, it is not clear whether the snapshots that have been extracted are captured before, during or after the peak of the event, dened as the time at which the local Hölder exponent reaches its minimum.

Link with the energy transfers

The local energy transfer D I stays positive in the averaged event, even though strong negative energy transfers are observed in individual events. This means that while there are large uctuations, on average, extreme singular events transfer energy to the smaller scale. Furthermore, we observe that, at the scale = 1.8η, the average of the dissipative term D ν in the typical event is slightly stronger than the average of the energy transfer term D I . Using the discussion from Chapter 2, this means that events used in this average are observed approximately at the scale at which they get regularized by viscosity. Note that as the study has been performed very close to the Kolmogorov scale, the Duchon Robert term D I is likely still comparable to the dissipative term D ν at = η. This implies a contribution to the energy budget which would not have been captured if the cuto scale for the simulation was the Kolmogorov scale.

Another observation that can be made on the averaged event is that the structures of D I and D ν are not on the same axis. From the point of view of the streamlines, the peak of energy transfer happens before the dissipation. This is not the eect of the averaging, as such structure is also observed in individual events (see gure 4.22).

Temporal evolution of an event.

The purpose of using a time resolved dataset is to try and uncover the generation mechanism for extreme events. In the time period covered by the time resolved part of the simulations, we found 11 events for which the minimum of local Hölder exponent is kept under the threshold of h = 1.3 for 10 snapshots or more. This corresponds to at least 3 eddy turnover times t η . This very strict threshold keeps the number of events to analyze low. Indeed, analyzing time series of events is a lot more time consuming than what was done in the previous section.

The results support the hypothesis evoked previously: the extreme events of local Hölder exponent appear related to vortex interaction. Indeed, two or more vortices are observed for most of the events detected. One such event is illustrated at dierent times in gures 4.27 and 4.28. For this event, the structures of Hölder, energy transfer and energy dissipation start to appears in gure 4.27b, when the vortices are approaching each other. We recognize the structure of the typical event in gure 4.27d, as the vortex observed is deformed by another vortex coming from z > 0. This structure becomes more chaotic when the vortices further interact in gure 4.27e. For this particular event, the structure of strong energy dissipation D ν does not follow the structures of low Hölder exponents h and high energy transfer D I in the second part of the interaction visualized in gure 4.28. The structures shown in gure 4.28 dier from the typical event obtained by averaging in gure 4.24. This might be the origin of the complex structures identied in the uncorrelated dataset. It is worth noting that the local Hölder exponents becomes greater than 1.3 before the reconnection is completed.

A timeline of the evolution of the scalars analyzed is given in gure 4.26.

This event is particularly irregular, as the local Hölder exponent drops to almost h = 1.2. The rarity of the event can be visualized by how low this value is on gure 4.15.

From the time evolution of the event, the most singular moment happens during the interaction, as the vortices are close to each other, at the time corresponding to gure 4.27f. At this time, the vorticity and the energy dissipation reach their maximum, and the Hölder exponents reach their minimum. Interestingly, we observe that the energy transfer D I reaches its maximum slightly before the other scalars. This behavior is observed for multiple events. This coincides with the previous observation that regions of high energy transfer are located earlier on the streamlines of velocity. While this could be a hint toward the mechanism behind the formation of quasi singularities, we currently lack a physical explanation for Chapter 5

A PIC simulation method using the Clebsch decomposition

The method exposed in the previous chapters allows us to localize potential singularities of NavierStokes in turbulent velocity elds. Because of the relatively low Reynolds number at which the method has been used, we are unable to follow those structures very far in scale. The objective of this chapter is to try and develop a new numerical method for NavierStokes simulations based on a Particles in Cells (PIC) method. This would allow us to use much smaller viscosity and potentially follow singular events much farther in scale. The general construction of this chapter is widely inspired from the thesis of Jean-Philippe Laval [START_REF] Laval | Développement d'un nouveau modèle dynamique pour la turbulence[END_REF], who developed a PIC method for NavierStokes. The original method relies on particles of vorticity advected by the large scale velocity in order to model the small scales of the ow. In this chapter, we use the Clebsch decomposition of the velocity eld to exhibit scalars eld that can fulll the same role as vorticity did. We then lay the basis for the simulation model using those scalars for the PIC method.

Motivations

In addition to the considerations for the simulations of singularities, the trigger for the work performed in this chapter comes from three publications, respectively by Yoshida and Morrison [START_REF] Yoshida | Epi-two-dimensional uid ow: a new topological paradigm for dimensionality[END_REF], Cartes et al. [12] and Ohkitani and Constantin [START_REF] Ohkitani | Eulerian-lagrangian analysis of navier-stokes turbulence[END_REF].

Among those, the rst one is dedicated to an analysis of the topology of the ow, while the two others develop simulation methods for NavierStokes. The common point of these papers is the usage of a Clebsch parameterization of the velocity eld, also known as the Weber-Clebsch representation. This method decomposes the initial vector eld into several scalar elds, for which the evolution equations can be derived. Using the reasoning reproduced in section 5.4.1, these papers show that the evolution equations for these scalar elds are similar to the equation of diusion of a passive scalar. This makes this scalar elds good candidate to the use of a PIC method: we expect to be able to model the small scales of the ow by particles carrying the information of the scalar elds passively advected by the large scale ow. This bears great similarity with the work performed by

Laval [START_REF] Laval | Développement d'un nouveau modèle dynamique pour la turbulence[END_REF]. In his thesis, Laval models 2D turbulence using a PIC method using the equation for vorticity. However, the approach used can not be generalized to the 3D equations of vorticity. The advantage of a method using the Clebsch parameterization over this previous work is that it would ba applicable to the 3D NavierStokes equation.

Decomposition between large scales and small scales

As a prerequisite to develop the new PIC method, we recall the NavierStokes equations in vector form:

∂ t u + (u • ∇)u = - 1 ρ ∇p + ν∆u + f ∇ • u = 0 (5.1) (5.2)
As expressed before, the general idea behind the PIC method is to decompose the velocity eld into large scales and small scales. The large scales are to be simulated using a classical numerical method, while the small scales are to be approximated by particles advected by the large scales. Formally, the decomposition is as follow:

U (x) = ū(x) = G(x -x )u(x )dx , u (x) = u(x) -U (x).
(5.3)

(5.4)

Applying the lter on equation (5.1) yields:

∂ t U + (u • ∇)u = - 1 ρ ∇P + ν∆U + F (5.5)
In a way not dissimilar to Large Eddy Simulations, this equation will be integrated while approximating the term (u • ∇)u thanks to the small scales. For this purpose, we get the small scale equation by subtracting equation (5.5) from equation (5.1). This leads to:

∂ t u + (u • ∇)u -(u • ∇)u = - 1 ρ ∇p + ν∆u + f (5.6)
At this point we need to make some hypothesis in order to simplify the equations. Those hypothesis are based on two little parameters. The rst one is based on the ratio between the large scales L and the small turbulent structures l: = l/L. The other one is based on the ratio between the characteristic time associated to the turbulent scales τ and the characteristic time associated to the large scales T (estimated as T = 1/∇U ) : η = τ /T . Using the results from Dubrulle

Nazarenko [START_REF] Dubrulle | Interaction of turbulence and large-scale vortices in incompressible 2d uids[END_REF], we take the hypothesis 1/η 1.

Applying dimensional analysis to the non linear terms and taking (U • ∇)U as a reference, we have in orders of magnitude:

(U • ∇)u ∼ 1/( η) (U • ∇)U (u • ∇)U ∼ /η (U • ∇)U (u • ∇)u ∼ 1/η 2 (U • ∇)U
Keeping only the dominant terms, equations (5.5) and (5.6) reduce to:

∂ t U + (U • ∇)U = - 1 ρ ∇P + ν∆U + F f orcing + F small→large ∂ t u + (U • ∇)u = - 1 ρ ∇p + ν∆u + F large→small (5.7)
(5.8)

with:

F small→large = -(U • ∇)u F large→small = (U • ∇)U -(U • ∇)U + (U • ∇)u .
(5.9)

(5.10)

The approximations realized here have already been veried in [START_REF] Laval | Développement d'un nouveau modèle dynamique pour la turbulence[END_REF].

On the other hand, applying the decomposition between large and small scales to the incompressibility equation (5.2) yields:

∇ • U = 0 and ∇ • u = 0.
(5.11)

(5.12)

Algorithm for the large scales

For this model, the large scales can be simulated using a conventional algorithm, similar to what is used for LES simulations. Our choice in this case is to use a spectral code, as it has the merit of simplicity.

In order to derive the equations for the spectral method, we take the Fourier transform and solenoidal projection of (5.7). Assuming the system is forced directly in the spectral space, this results in the following:

∂ t U (k, t) = -νk 2 U (k, t) + G(k, t) + F f orcing (k, t) (5.13)
where G(U ) gathers the non linear terms:

G(k, t) = F T (P sol (-(U • ∇)U + F small→large )) (5.14)
where F T designate the Fourier transform and P sol the solenoidal projection.

Putting the forcing aside and introducing U (t) = e νk 2 t U (t), we get:

∂ t U (k, t) = e νk 2 t G(k, t). (5.15)
This equation can be integrated numerically. One possible choice is to use a Runge-Kutta algorithm. For a time step of ∆t, this would give:

U (k, t + ∆t) -U (k, t) = ∆te νk 2 (t+∆t/2) G(k, t + ∆t/2).
(5.16)

Going back to the original variables, the integration scheme is the following:

U (k, t + ∆t) = e -νk 2 ∆t U (k, t) + ∆te -νk 2 ∆t/2 G(k, t + ∆t/2) (5.17)
where the computation of G(k, t + ∆/2) will require the use of:

U (k, t + ∆t/2) = e -νk 2 ∆t/2 U (k, t) + ∆t 2 G(k, t) .
(5.18)

As a result, the core of the model consists in evaluating the non linear term G with a good balance of accuracy and eciency. This prominently factor the retro-action from the small scales to the large scales F small→large .

The forcing can be added every time step simply as follow:

U f orced (k, t + ∆t) = U (k, t + ∆t) + ∆tF f orcing (k, t).

(5.19)

Equations for small scales

In this section, we derive the equations governing the small scales that will be implemented in a PIC model. This will eventually lead to the computation of the small scale contribution to the non linear term G.

Clebsch decomposition

Deriving the equations with viscosity massively increases the complexity of the computation. This work has been done, in particular in [12], but the resulting equations are less suitable to a PIC method than in the inviscid case. For the time being, we consider the inviscid equation without the large scale contribution:

∂ t u + (U • ∇)u = - 1 ρ ∇p .
(5.20)

Making the Particles of our Particles in Cells model not dissipative also participates in our overall objective of trying to observe singular behaviors. Indeed, this means that the small scales quantities will follow Euler's equations, and potentially transport energy to arbitrarily small scales.

The Clebsch decomposition is introduced by: u = -∇φ + p α=1 λ α ∇µ α .

(5.21)

In this equation, φ is a scalar eld working exactly as the potential in the Helmholtz decomposition. The specicity of the Clebsch decomposition comes from the scalar elds λ α and µ α for α = 1..p. As we are working in 3 dimensions, a result by Yoshida [START_REF] Yoshida | Clebsch parameterization: basic properties and remarks on its applications[END_REF] guarantees the existence of such scalar elds for p 2. There is however no unicity of the decomposition.

Injecting the equation (5.21) into (5.20) gives:

p α Dλ α Dt ∇µ α - Dµ α Dt ∇λ α -∇ Dφ Dt - p α Dµ α Dt λ α -u • (∇U ) = - 1 ρ ∇p (5.22)
where we use the material derivative D Dt = ∂ t + (U • ∇), computed using the large scale velocity. In this equation, the term colored in red is of the same order as (u • ∇)U which has already been neglected in the previous section, so we neglect it again.

For this equation to be satised, a sucient but not necessary condition is:

Dλ α Dt = 0 Dµ α Dt = 0 Dφ Dt = p ρ (5.23)
(5.24)

(5.25)

Combined with the non unicity of the scalar elds λ α and µ α in the Clebsch decomposition, we can aord to impose these conditions. The equations (5.23) and

(5.24) are particularly interesting as it means the scalar elds can be portrayed as "passively" advected by the large scale velocity. This description is inaccurate as the model we develop here does include a retro-action on the large scale velocity, as described in the section 5.5.5. In this aspect, it is a situation strikingly similar to how the small scale vorticity is "passively" advected by the large scale in Laval's thesis [START_REF] Laval | Développement d'un nouveau modèle dynamique pour la turbulence[END_REF]. This similarity is what lead us to try and develop a PIC based method, where the scalar elds λ α and µ α are decomposed into particles advected by the large scale velocity.

Supposing that the scalar elds λ α and µ α are known, we can recover the potential eld φ using the incompressibility condition. Indeed, applying the Clebsch decomposition from equation (5.21) to the equation (5.12) yields: ∆φ = p α ∇ • (λ α ∇µ α ). (5.26) This means we can avoid integrating the equation (5.25), and thus that we don't have to keep track of the small scale pressure.

Filter function and Gabor transform

In order to progress further, the ltering needs to be formalized. The lter is dened as:

U (X, t) = ū(x, t) = f 2 ( (x -x )) u(x )dx (5.27)
where is a small parameter dierentiating the "slow" spatial variable X = x/ and the "fast" variable x. The small parameter is such that 1

.

The function f is an odd function with a fast decay toward innity and a norm unity, i.e.

f 2 (x)dx = 1. This ensures an intuitive behavior of the ltering when applied to constant functions. The choice of a f 2 instead of f does not change the properties of the problem and will simplify later computations.

Additionally, we introduce the Gabor transform using the same function f :

GT {u }(x, k, t) = û (x, k, t) = f ( (x -x
)) e ik(x-x ) u (x , t)dx .

(5.28)

This operation bears striking similarities with continuous wavelet transforms, as it conveys information about both space and frequency. It is the principal factor by which the small scale scalar elds will be modelled by particles evolving in the space-frequency domain.

Similarly to the continuous wavelet transforms, this transformation can be inverted as follow:

u (x, t) = 1 (2π) 3 f (0) û (x, k, t)dk.
(5.29)

All elds in physical space are real valued, which leads to:

û (x, -k, t) = û (x, k, t) * λα (x, -k, t) = λα (x, k, t) * μα (x, -k, t) = μα (x, k, t) * (5.30)
where the star denotes the complex conjugate.

Several properties of the Gabor transform will be required for later computations. Those properties are also demonstrated in [START_REF] Laval | Développement d'un nouveau modèle dynamique pour la turbulence[END_REF] and these demonstrations are reported here for the sake of convenience. By using integration by part, one can show that the Gabor transform and spatial derivatives commute:

T G{∂ i u }(x, k, t) = ∂ i û (x, k, t).
(5.31)

The same term can also be expanded using the small parameter . Indeed, we have:

∂ i û (x, k, t) = ∂ i f ( (x -x )) e ik(x-x ) u (x , t)dx = ik i û (x, k, t) + f ( (x -x )) e ik(x-x ) u (x , t)dx = ik i û (x, k, t) + O( ) (5.32)

Main equations for the small scale velocity

The objective in this part is to derive the equations for the Gabor transformed scalar elds λα and μα which would then be used as basis for the PIC model.

For this purpose, we apply the Gabor transform to equation (5.23). In the following Einstein summation convention for repeated indices is used:

∂ t λα + GT {U j ∂ j λ α } = 0.
(5.33)

We need to compute the Gabor transform of a product of a small scale quantity and a large scale quantity:

GT {U j ∂ j λ α }(x, t) = f ( (x -x
)) e ik(x-x ) U j (x , t)∂ j λ α (x , t)dx .

(5.34)

As the large scale quantity U is slowly varying compared to f , we perform a Taylor expansion of U around x :

U j (x , t) = U j (x, t) + (x -x ) • ∇U j (x, t) + O((x -x) 2 ).
GT {U j ∂ j λ α }(x, t) ≈U j (x, t)GT {∂ j λ α }(x, t) + ∂ l U j (x, t) f ( (x -x )) e ik(x-x ) (x l -x l )∂ j λ α (x , t)dx ≈U j (x, t)GT {∂ j λ α }(x, t) + i∂ l U j (x, t) ∂ ∂k l GT {∂ j λ α }.
(5.36)

Finally, using the propriety from equation (5.32), we have:

0 = ∂ t λα + GT {U j ∂ j λ α } ≈ ∂ t λα + U j ∂ j λα + i∂ l U j ∂ ∂k l (ik j λα ) ≈ ∂ t λα + U j ∂ j λα -∂ l U j k j ∂ ∂k l λα - ¨¨¨λ α ∂ j U j (5.37)
The term ∂ j U j is removed as a consequence of incompressibility from equation (5.11).

The exact same process can be applied to the equation (5.24). As a result, we have the main equations for the Gabor transformed scalar elds λα and μα :

D t λα (x, k, t) = 0, D t μα (x, k, t) = 0, (5.38)
(5.39) with :

D t = ∂ t + ẋ • ∇ + k • ∇ k , ẋ = ∇ k Φ, k = -∇Φ, Φ = U • k.
(5.40)

(5.41)

(5.42)

(5.43)

Implementation of a PIC method

The dierential operator D t from equations (5.38) and (5.39) can be interpreted as a material derivative in the space-frequency domain with a well chosen "velocity" directly dened from the large scale velocity U . This way, the Gabor transformed scalar elds λα and μα are interpreted as scalar transported in the space-frequency domain. The idea of the PIC model is to decompose the small scale scalars into particles with given position and wave vector that are advected independently using equation (5.40). The retro-action of the small scale on the large scale motion would then be a combination of the action of all individual particles, with the estimation being more precise the more particles there are.

Such computing method have been rst developed for computing the dynamic of charged particles within a plasma and has been used eectively ever since [8]. Its use has been expanded to less straightforward application with no actual particle involved, for example by Nazarenko, Zabusky and Sheidegger in order to simulate the interaction between wave packet of sound and vortices [START_REF] Nazarenko | Nonlinear soundvortex interactions in an inviscid isentropic uid: A two-uid model[END_REF]. It has also been used by Laval to simulate the transport of vorticity [START_REF] Laval | Développement d'un nouveau modèle dynamique pour la turbulence[END_REF].

The expected advantage of such method over a direct numerical simulation is that an individual particle might represent much smaller scale for a meager computation cost compared to a renement of the resolution.

In order to implement the PIC model, we will rst discuss the denition of the particles and how they are related to the scalar elds λα and μα . We will then proceed to detail how those particles implement the dierent contributions of small scales on the large scales of the ow. This means implementing the forcing F large→small from equation (5.10) that large scales exert on small scales in the form of creation of new particles. We will then discuss the implementation of a

k 1 λα,1 μ α,1 k 2 λα,2 μ α,2 k 3 λα,3 μ α,3 x 1 x 2 U(x 1,1 ) U(x 2,1 ) U(x 3,1 ) U(x 3,2 ) U(x 2,2 ) U(x 1,2 ) U(x 1,3 ) U(x 2,3 ) U(x 3,3 )
x 3 ∆x numerical dissipation based on the removal of particles from the computations. The description of the implementation of the advection of particles will follow.

Finally, we will explain how to derive the retro-action of the small scales on the large scales F small→large using the particles.

Denition of the particles

Let us suppose that the small scale scalar elds λα et μα are decomposed into N p particles characterised by:

          
x n the posistion of the particle, k n its corresponding wave vector (3D), λα,n the complex amplitudes for the λ elds, μα,n the complex amplitudes for the µ elds.

The decomposition has to be such that the scalar elds can be reconstructed as follow:

λα (x, k, t) = Np n=1 λα,n S x (x -x n (t))S k (k -k n (t)).
(5.44)

where S x et S k are the interpolation functions for the PIC method.

By simplicity, we chose the interpolation function in the wave vector space to be the Dirac function: S k (k) = δ(k). This zero order function will spare us complications when dealing with particle interactions. For the interpolation function in space, we choose a simple rst order interpolator:

S x (x) = S x (x) × S y (y) × S z (z), (5.45)
where

S x (x) = dh -|x|/dh if |x| < dh, 0 otherwise.
(5.46)

and likewise for S y and S x . We choose the parameter dh equal to the grid step for the large scales dh = ∆x.

The gure 5.1 gives a representation of particles moving independently from the main simulation grid used for the large scales. In this representation, the eld λα at the point x 2,2 can be reconstructed from the complex amplitudes carried by the particles numbered 1,2 and 3. With this choice of interpolation function S x , the reconstruction of a small scale eld at a given node only depends on particles directly around the node. A higher order of interpolation might make the result more accurate, but the computational cost will be heavier.

One can observe that the equations (5.38) and (5.39) are symmetric under the transformation k → -k.

For our PIC model, it means that for every particle simulated, one can deduce the behaviour of a particle with the same position but opposite wave vector. In practice, it means we can simulate twice as many particles for the same cost. Let us assume we add N p more particles to the system. We will distinguish pairs of particle with same position but opposite wave vector with a plus or minus index.

The equation (5.44) then becomes:

λα (x, k, t) = Np n=1 λα,n+ δ(k n (t)) + λα,n-δ(-k n (t)) S x (x -x n (t)).
(5.47)

By taking the inverse Gabor transform of this equation using the formula (5.29), this results in:

λ α (x, t) = 1 (2π) 3 f (0) Np n=1 λα,n+ + λα,n-S x (x -x n (t)).
(5.48) Because the resulting scalar elds are real valued, this imposes λn,α+ = λ * n,α-.In the following, we drop the + index when there is no ambiguity. Adopting the notation [z] and [z] for respectively the real part and the imaginary part of a complex number z, the resulting equation is:

λ α (x, t) = 2 (2π) 3 f (0) Np n=1 λα,n S x (x -x n (t)).
(5.49)

Of course, the exact same reasonning can be applied for µ, leading to:

µ α (x, t) = 2 (2π) 3 f (0) Np n=1 [μ α,n ] S x (x -x n (t)) (5.50)
When a spatial derivative is used, the property of the Gabor transform with respect to derivation (5.32) gives:

∂ j µ α (x, t) ≈ -2 (2π) 3 f (0) Np n=1 k j [μ α,n ] S x (x -x n (t)).
(5.51)

Forcing from the large scales: creation of particles

Similarly to the PIC method developed in the thesis of Laval, the eect of large scales on small scales is not modeled by interacting with the existing particles but by creating new particles. At this point, the equation for the small scales is as follows:

D t u = - 1 ρ ∇p + F large→small .
(5.52)

The pressure term is already included in the dynamic of the particles. In order to model the contribution of the large scale to the small scales, we add at each time step particles corresponding to the velocity eld u f = ∆t × F large→small , where ∆t is the time step for the large scale simulation dened in section 5.4. This is the opportunity to decide the number p of scalar elds used in the Clebsch decomposition, which has been left undetermined until now. We chose the value p = 3, which is also the value taken by Cartes et al. [12]. As this corresponds to the number of spatial dimensions, we adopt the same convention for the indices in the following.

This means that we can take the same choice as [12] for the Clebsch decomposition of u f :

λ i (x) = u f i , µ i (x) = x i .
(5.53) (5.54) This particular choice is inconvenient when using periodic boundary conditions.

For this purpose, we set µ i (x) = µ i (x)x i and the complex amplitudes carried by the particle model will correspond to this new quantity µ . Note that even if the µ i are not periodic, their derivatives ∂ j µ i are, so there is no problem for the decomposition of u f periodic.

From this change of variable, it entails that:

μ i (x, k) = μi (x, k) -xi (x, k).
(5.55)

This leads to:

∂ j μ i (x, k) = ∂ j μi (x, k) -δ ij 1(x, k) = ∂ j μi (x, k) -δ ij f (k) (5.56) (5.57)
where 1 is the Gabor transform of the function unity, which is also the Fourier transform of the lter function f (k).

The equation (5.39) becomes:

D t μ i (x, k, t) = -U i f (k).
(5.58)

The elds λ and µ corresponding to this added velocity u f being determined, we needs to decompose these elds into particles to add to the system. A simple choice is to generate N 3 particles, indexed by n, with one particle at each node of the grid used for the simulation of the large scales. With the initial complex amplitude corresponding to the µ elds being set at zero, we still need to determine the wave vectors k n and the complex amplitudes λ i,n .

Applying equation (5.49) gives:

u f i (x) = 2 (2π) 3 f (0) N 3 n=1 λα,n S x (x -x n ).
(5.59) By evaluating this formula at each grid point, the support of the interpolation function S x is such that we have:

λi,n = (2π) 3 f (0) 2 u f i (x n ).
(5.60)

In order to obtain constraints on the wave vector k n and the imaginary part of the complex amplitudes, one must use the derivatives of the velocity eld u f . Note that if we use a spectral method to simulate the large scales, those derivatives are exact quantities. By taking two components of the rotational ω f = ∇ × u f as well as the incompressibility condition, we obtain the following system:

     0 λ3,n -λ2,n -λ3,n 0 λ1,n λ1,n λ2,n λ3,n         k1 k2 k3    =    -(2π) 3 f (0) 2 ω 1 -(2π) 3 f (0) 2 ω 2 0    (5.61)
Further derivatives might give more constraints, at the cost of increased complexity. We follow the choice made by Laval in similar circumstances in the model described in [START_REF] Laval | Développement d'un nouveau modèle dynamique pour la turbulence[END_REF], and impose: λi,n = λi,n (+ )

(5.62)
where the small quantity is added when necessary to prevent the matrix in (5.61) from being non inversible.

To summarize, the N 3 particles created at each time step at all node of the large scale grid have the following initial parameters:

                
x n the position of the grid point, k n the wave vector solution of (5.61), λi,n = (2π) 3 f (0) 2 u f i real part of the amplitude for λ, λi,n = λi,n (+ ) imaginary part of the amplitude for λ, μ i,n = 0 initial complex amplitude for µ set to zero,

(5.63)

The generation of new particles in this manner causes several complications for the model. The rst one is caused by the properties of the Clebsch decomposition.

Indeed, there is no unicity for the elds λ and µ. The concrete consequence is that the λ and µ elds decomposed at time T cannot be combined with the elds at time T + ∆T . In other words, at each time step the Clebsch decomposition of the velocity eld u f is performed using a dierent gauge condition.

This diculty can be overcome by keeping track of the time at which the particles have been created. Assuming that the elds λ T i (t) and µ T i (t) are recovered from the particles created at time T , the small scale velocity at time t is computed as such:

u (x, t) = T <t u T (x, t) = T <t -∇φ T (x, t) + λ T α (x, t)∇µ T α (x, t) .
(5.64)

The computation is then developed as follow. First, the incompressibility con-dition on the Clebsch decomposition from equation (5.26) is used to get rid of the potential elds. We then change from the µ elds to the modied µ :

u j (x, t) = δ jl -∂ j ∆ -1 ∂ l T <t λ T α (x, t)∂ l µ T α (x, t) = δ jl -∂ j ∆ -1 ∂ l T <t λ T α (x, t) ∂ l µ T α (x, t) + δ αl (5.65)
The nal expression is obtained by injecting the equations (5.49) and (5.51):

u j (x, t) = 4 (2π) 6 f (0) 2 δ jl -∂ j ∆ -1 ∂ l T <t N T n=1 λT α,n S x (x -x n (t)) δ jl - N T m=1 k l μ T α,m S x (x -x m (t)) (5.66)
where the sums over n and m are performed over the particles created at time T .

In this expression, the dierential operator preceding the sum can be computed using a spectral method over the large scales. It is worth noting that this operation's computational complexity is quadratic in the number of particles created at time T , i.e. in N 3 . While the small support of the interpolation function keeps the actual number of operation moderate, we will design the algorithm in a way that minimizes the use of this formula.

Dissipation

The other problem caused by the generation of particles is of practical order.

Indeed, running such model would imply adding N 3 particles to the system at each time step. Following through this process would quickly increase the cost of computations, such that the computational advantage over a more classical method would be lost. We list here several ideas that could be used to keep the number of particles under control.

The rst idea is to use the removal of particles to implement the numerical dissipation. The most natural approach would then be to use the equivalent of a spectral lter, removing all particles with a wave vector too large k > k max . This is coherent with the method used by Laval in his thesis [START_REF] Laval | Développement d'un nouveau modèle dynamique pour la turbulence[END_REF].

Our situation present one major dierence with the case discussed in [START_REF] Laval | Développement d'un nouveau modèle dynamique pour la turbulence[END_REF].

This distinction is that we use the particle model to model the small scale velocity instead of the small scale vorticity. This remove the possibility to create particles only in the representative points, as there is no reason to not simulate particles based on a threshold on velocity. We could use a threshold based on vorticity, at the cost of computing the rotational every time particles are created.

One option in common with [START_REF] Laval | Développement d'un nouveau modèle dynamique pour la turbulence[END_REF] is the possibility to "recycle" the particles after a certain time. The idea is that once the number of particles has reached a certain threshold N max , we can use equation (5.66) to reconstruct the small scale velocity u . This reconstructed velocity can then be transformed back into particles using the same process as described in part 5.5.2. In our case, this method can be applied to only "recycle" part of the particles by delimiting a range of time.

In other words, we could choose to "recycle" only the particles from time T 1 to T 2 . This would ensure that the number of particle do not fall so low that the PIC reconstruction loses in precision.

Recycling the particles also allows to limit the risk of a particle falling into a local reverse cascade, with the norm of the wave vector k reaching a value too low.

This would cause the formalism to break as the small scale velocity would have a slowly varying component. While this event is unlikely in 3D turbulence, it is not forbidden locally.

All those methods require testing to verify the eect on the accuracy of the simulation and on the computational cost. As was discussed in the introduction to this chapter, the implementation of the algorithm has not been completed yet, hence nether are the tests.

Advection of particles

The dynamic of the particles is governed by equations (5.38) and (5.58). Expressed from the point of view of the particles parameters, these translate to:

∂x n ∂t (t) = U (x n (t)), ∂k n ∂t (t) = -(k n (t) • ∇)U (x n (t)), ∂ μ i,n ∂t (t) = -U i f (k n (t)).
(5.67)

(5.68)

(5.69)

Incidentally, the complex amplitude associated to the λ eld is constant:

λi,n = cst.

(5.70)

In order to implement these equations, we need to know the large scale velocity eld U and its derivative at a point which is not on the grid. For this, we need to use the already dened interpolation function S x from equation (5.45).

Formally, we consider the 8 closest nodes to the particle n, noted x n,β for β = 1 . . 8. Then the large scale velocity at the location x is obtained by:

U (x, t) = 8 β=1 U (x n,β , t)S x (x n,β -x n ).
(5.71) Going back to the schematic representation of gure 5.1, one can visualize how the velocity at the position of a particle is computed as the weighted average of the velocity at the neighboring grid nodes. Like for the reconstruction, it is possible to choose a higher order for the interpolation function S x in order to get more accurate results but this would be at higher computational cost.

The same can be done for the dierent spatial derivatives of U , provided that these are computed exactly using the spectral method on the large scale grid.

The equations (5.67) to (5.69) can then be integrated using the second order stable Runge-Kutta method. In order to compute the new parameters of the particles at a time step ∆t/2 required by the large scale scheme from equation (5.18) we need the intermediate variables:

x n = x n (t) + ∆t 4 U (x n (t)), k n = k n (t) - ∆t 4 (k n (t) • ∇)U (x n (t)).
(5.72)

(5.73)

The time stepping over the duration ∆t/2 is then given by:

x n t + ∆t 2 = x n (t) + ∆t 2 U (x n ), k n t + ∆t 2 = k n (t) - ∆t 2 (k n • ∇)U (x n ), μ i,n t + ∆t 2 = μ i,n (t) - ∆t 2 U i (x n )f (k n ).
(5.74)

(5.75)

(5.76)

Note that the currently chosen numerical scheme means performing operations on the particles four times by time step in total.

Retro-action on large scales

Having determined the mechanic of creation/removal of particles and their dynamic, we can use them to determine the contribution of the small scale velocity to the large scales. The term involving the small scale velocity that need to be determined is U i ∂ i u j , which is present for both the forcing on small scales F large→small (equation (5.10)) and the retro-action on large scales F small→large (equation (5.9)).

To simplify this expression, we use the Taylor development of the large scale velocity U :

U i ∂ i u j (x, t) = f 2 ( (x -x )) U i (x , t)∂ i u j (x , t)dx ≈ f 2 ( (x -x )) [U i (x, t) + (x -x ) • ∇U i (x , t)] ∂ i u j (x , t)dx ≈ U i (x, t)∂ i u j (x, t) + f 2 ( (x -x )) (x -x ) • ∇U i (x , t)∂ i u j (x , t)dx = U i (x, t)∂ i u j (x, t) + O( ).
In this last development, only the rst order is kept. This is because the next term is of comparable magnitude with terms already neglected when deriving the equations in section 5.2. As a result, we have:

F small→large ≈ -(U • ∇)u .
(5.77)

The nal step consists in determining the ltered small scale velocity u (x, t)

using the particle model. One possible option would be to compute the small scale velocity u (x, t) using the formula from equation (5.66) then apply the lter.

However, it has already been mentioned that this operation is computationally costly as it is quadratic with respect to the number of grid points N 3 . A cheaper alternative is to use the following:

u j (x, t) = T <t u T j (x, t) = δ jl -∂ j ∆ -1 ∂ l T <t λ T α ∂ l µ T α (x, t).
(5.78)

In order to compute the lter of a product of two small scale quantities, we make the following remark:

λT α (x, k, t)∂ l μT α (x, -k, t) dk = f ( (x -x ))e ik(x-x ) λ T α (x , t)dx f ( (x -x ))e ik(x-x ) ∂ l µ T α (x , t)dx dk = f ( (x -x ))f ( (x -x ))λ T α (x , t)∂ l µ T α (x , t) e ik(x -x ) dk dx dx =(2π) 3 λ T α ∂ l µ T α (x, t).
(5.79)

This leads to:

λ T α ∂ l µ T α (x, t) = 1 (2π) 3 λT α (x, k, t)∂ l μT α (x, -k, t) dk = 1 (2π) 3 λT α (x, k, t)(∂ l μ T α (x, -k, t) + δ lα f (-k)) dk ≈ 1 (2π) 3 λT α (x, k, t)(-ik l μ T α (x, -k, t) + δ lα f (-k)) dk ≈ 1 (2π) 3 λT l (x, k, t)f (-k) dk + 1 (2π) 3 k l λT α (x, k, t)μ T α (x, -k, t) dk.
(5.80)

The rst term can be computed directly by applying the relation (5.44). In the following computation, the rst sum needs to include the particles with opposite wave vectors introduced in part 5.5.1. The lter f being real valued, we use the properties f (k) = f (-k) * and λT l,n-= λT * l,n+ to simplify the expression.

λT l (x, k, t)f (-k) dk = N T p n=1 λT l,n+ δ(k -k T n (t)) + λT l,n-δ(k + k T n (t)) f (-k)S x (x -x T n (t)) dk = 2 N T p n=1 λT l,n f (-k T n ) S x (x -x T n (t)).
(5.81)

The second term of equation ( 5.80) can be computed using the same method.

k l λT α (x, k, t)μ T α (x, -k, t) dk = k l     N T p n=1 λT l,n+ δ(k -k T n (t)) + λT l,n-δ(k + k T n (t)) S x (x -x T n (t))     N T p m=1 μ T l,m+ δ(-k -k T m (t)) + μ T l,m-δ(-k + k T m (t)) S x (x -x T m (t))     dk.
Using the fact that dierent particles are very unlikely to have the exact same wave vector, the choice of a Dirac for the interpolation function in the spectral space in part 5.5.1 allows us to simplify the computation:

k l λT α (x, k, t)μ T α (x, -k, t) dk = 2 N T p n=1 k T l,n λT α,n μ T * α,n S x (x -x T n (t)) 2 .
(5.82)

Finally, injecting the results of (5.81) and (5.82) in equation ( 5.80) and grouping the sums, the ltered small scale velocity is recovered as follow:

u j (x, t) = 2 (2π) 3 δ jl -∂ j ∆ -1 ∂ l Np n=1 λl,n f (-k n ) S x (x -x n (t)) +k l,n λα,n μ * α,n S x (x -x n (t)) 2 .
(5.83)

The leading dierential operator of this expression can be computed using a spectral method. Remarkably, the computation cost of this term is linear in the number of particles.

Constraints on the lter

The knowledge of the lter is required in order to complete the model. However, the choice of this lter is not left free. This constraint is apparent in the computation of the lter of a product of small scale variables. For the sake of argument, let us take the ltered product λ T i µ T i . A common rectangle integration over all the grid points x i would give:

λ T i µ T i (x, t) = i λ T i (x i )µ T i (x i )f 2 (x -x i )(∆x) 3 .
(5.84)

Estimating the same quantity using the particle model yields:

λ T i µ T i (x, t) = 2 (2π) 3 Np n=1 λT i,n μT i,n S x (x -x T n ) 2 .
(5.85)

As those two relations must be true for any two small scale elds, this imposed that f 2 and S 2

x are proportional to one another:

f 2 = C × S 2
x .

(5.86)

Additionally, the lter must be normalized such that f 2 (x)dx = 1. As a result, we have:

f 2 = f 2 (0)S 2 x .
(5.87)

where

f (0) = 1 S 2 x (x)dx
.

(5.88)

For our choice of S x , this results in:

f (0) = 3 2dh 3/2 . (5.89)
The Fourier transform of both the function and its square are required. Those are computed as follow:

f (k) = f (0)dh 3 sinc 2 k x dh 2 sinc 2 k y dh 2 sinc 2 k z dh 2 , f 2 (k) = 6 3 sinc(k x dh) -1 (k x dh) 2 sinc(k y dh) -1 (k y dh) 2 sinc(k z dh) -1 (k z dh) 2 .
(5.90)

(5.91)

As a remark, we have:

lim x→0 sinc(x) -1 x 2 = - 1 6 .
(5.92) so there is no divergence around the origin for f 2 .

Accessing to small scale statistics: the energy spectrum

Assuming the algorithm described in this chapter is implemented, we would need to have some access to statistical values pertaining to the small scales. At the very least, this would be required in order to validate the model. At best, we could make use of these statistics to study some small scale behaviors with only a fraction of the computational cost of a large DNS.

The very rst quantity that needs to be checked is the energy spectrum. Again, assuming the simulation runs smoothly and the energy spectrum of the large scale behave as expected, showing an inertial range, we would then focus our interest on the energy spectrum of the small scales. As the equations for the dynamic of the particles are inviscid, we would expect the inertial range to extend from the large scales to the cuto k max at which we remove the particles.

The similar work by Laval [START_REF] Laval | Dynamical modeling of sub-grid scales in 2d turbulence[END_REF] used the energy carried by individual particles in order to reconstruct the complete energy spectra. This method is not strictly applicable the current model because of the non linearity of the Clebsch decomposition: the local energy is not a simple sum of the "energy" carried by neighboring particles. Nevertheless, we can still use this approach as an estimation. Starting from the expression of the velocity combined from the particles created at time T :

u T i = (δ ij -∂ i ∆ -1 ∂ j )(λ T j + λ T α ∂ j µ T α ).
(5.93)

The λ and µ terms can be recovered using the particles by:

λ T i = 2 (2π) 3 f (0) n λT j,n S x (x -x T n ), ∂ j µ T α ≈ - 2 (2π) 3 f (0) n k T j,n μ T α,n S x (x -x T n ).
(5.94)

(5.95)

If we ignore the contributions from the other particles, the velocity at point x n associated to a single particle n is expressed as:

u T i,n ≈ 2 (2π) 3 f (0) λT i,n - 4 (2π) 6 f (0) 2 k T i,n λT α,n μ T α,n .
(5.96)

As a result, the energy associated to this particle would be:

E T n = 1 2 i u T i,n 2 .
(5.97)

In order to recover the energy spectrum for the small particles, the idea would then be to sum the contributions of all particles with a wave vector of same length.

More formally, the result would be:

E(k) = 1 dk T <t | k T n -k|<dk/2 E T n .
(5.98)

where the interval dk would be chosen large enough to have several particles per bin.

This method has the aw of circumventing the base principle of the PIC model, by neglecting all the terms that result in the combination of several particles. It is however the best approach that we have developed at the time of the redaction of this manuscript.

Current state and discussion

The implementation of this PIC method is currently at the debugging phase. As the project has been stalled for some time, it is unlikely to be someday concluded.

Should it continue, we would expect several roadblocks before completion.

From the point of view of the model, we expect that the process of creation of particles described in part 5.5.2 might require some revisions. The points with some degrees of freedom are the choice of the equations to constrain the wave vector k and the choice of imaginary part for the complex amplitudes λα,n . By analogy with the similar model developed in Laval thesis [START_REF] Laval | Développement d'un nouveau modèle dynamique pour la turbulence[END_REF], we expect these choices to be valid. However, they need to be properly validated.

While the exact method of validation is yet undetermined, one potential approach would be to run an a posteriori validation test against a well resolved DNS.

The principle of such validation would be to run the particle dynamic using the large scale velocity of the DNS as the large scales U , and verify if the particle model reproduces the same small scale velocity. Another idea would be to compare the retroaction of the small scales on the large scales computed in section 5.5.5 to a common LES model, for example the Smagorinsky model [START_REF] Smagorinsky | General circulation experiments with the primitive equations: I. the basic experiment[END_REF].

At this point, it is dicult to evaluate how this model would perform, let alone if it could outperform existing models. It however shows potential with regards to our original motivation of investigating the simulation of quasi singularities.

Indeed, the particle dynamics runs according to approximated Euler equations.

This means that in absence of numerical cuto, we would expect particles carrying energy to evolve such that their wave vector k would grow very large, corresponding to very strong gradient, without viscous dampening. By following the track of a group of particles past the frequency cuto, this could grant us access to quantities such as the rate of production of quasi singular events without having to use an expensive DNS followed by the complex computations presented in the rst part of this thesis.

Chapter 6

Conclusions and perspectives

The starting point for this thesis was the question of the regularity of the solutions to the NavierStokes equations. Our approach to this question was to look for irregular events in turbulent velocity elds which could be interpreted as potential singularities.

In a rst part, we developed a method based on multifractal analysis to estimate the local regularity of a turbulent velocity eld through the Hölder coecients. This method complements an existing criterion based on strong local energy transfers. That older method has already been used for the detection of potential singularities [START_REF] Kuzzay | New method for detecting singularities in experimental incompressible ows[END_REF][START_REF] Kuzzay | Investigations on the relevance of Onsager's conjecture in real incompressible turbulence[END_REF]. Using this preexisting concept as a reference, we validated in the inertial range the value of the method based on local Hölder exponents as a new criterion for the detection of irregular events. We indeed showed that the lower values of Hölder exponents did correspond to strong energy transfers, and thus to more irregular events. However, we also show by comparison with the VGT invariants that the criterion based on the local Hölder exponents is not limited to the same structures as the local energy transfer. This is conrmed by visualization of maps of these two quantities at the inertial scale.

We then described a way to use this tool in the dissipative range using a calibration in the inertial range. The results obtained in the dissipative range are in agreement with the previous study performed using experimental data [15].

The use of numerical simulation allowed us to produce statistical results without the diculties of experimental measurements. This comes with the drawback of a lower Reynolds number, which could potentially lead to less quasi singular events, or weaker events. Nevertheless, we obtained similar results for the topology of the streamlines of velocity around quasi singular events. Indeed, we observe the same shape of "roll vortices" previously identied by [15]. In addition, the amount of statistics available allowed us to dene a "typical event", with the same topology, 107 by averaging many very singular events. This shows the existence of a structure of local Hölder exponents, energy transfer and dissipation that is common to most irregular events. Furthermore, we exhibit a relationship between the apparition of potential singularities and vortex interactions. This is apparent from isolated snapshots of irregular events guring several vortices close to each other. This nding was conrmed by a detailed study of several events using a time resolved dataset. In particular, we found that the formation of a particular event coincides with an interaction between two vortices. Those ndings might be a step toward explaining the physical mechanism behind the formation of quasi singularities.

This in turn could potentially explain the phenomenon of anomalous dissipation.

In a second part, we tried to develop a new simulation scheme based on a PIC method. The principle of the method would be to model the small scale velocity as particles advected by the large scale. We show how this could be done using the Clebsch decomposition. We also derive a way to recover the eect of the small scale velocity on the large scales. The purpose of this new method of simulation in the context of this thesis would be to track potential singularities through the scales by the mean of the PIC model. The project is currently stalled before the new model could be validated.

6.1 Going forward

Universality of the local Hölder method

The method to compute local Hölder exponents developed in this thesis was meant to be later used on experimental data. This includes in particular the data from the study to which we compared our results in Chapter 4. However, we have shown that in order to apply the method to compute local Hölder exponents at the dissipative scale, one had to rst perform a calibration in the inertial scale.

This requirement to have a dataset with both an inertial range and a resolution close to the Kolmogorov scale is hard to fulll with experimental data. The solution to this problem has been mentioned in Chapter 3 and comes from the universal nature of turbulence. As such, the idea would be to determine the coecients c p corresponding to a new dataset without having to do the calibration in the inertial range.

First of all, we can verify that this universality exists. The paper by Geneste et al. [START_REF] Geneste | About universality and thermodynamics of turbulence[END_REF] shows that the structure functions verify an universal behavior. We can verify using our two simulations that this universal behavior is shared by the thresholds T p used to compute the active volumes, as expressed in equation (3.49). The gure 6.1a shows a good collapse of the curves of thresholds. For this gure, we chose to use the thresholds from run I as a reference, such that the corresponding β is β 1 = 1. In order to obtain the collapse, the β for run II is set to β 2 = 1.5. Moreover we observe that the coecients c p computed for the two datasets follow similar trends, as is shown in gure 6.1b. The work from here would be to compute the coecients c p for some other datasets in order to nd their dependency in Reynolds numbers. This would eventually allow us to apply the computation of local Hölder exponents to experimental data.

Improving the analysis

Obtaining the local Hölder exponents in experimental data, in particular with the dataset used by Debue [15] would be useful in determining whether the behaviors observed in Chapter 4 with simulations at relatively low Reynolds are robust.

We would be particularly interested in verifying that the typical structure of low regularity is the same at higher Reynolds number. It would also be interesting to see whether events with lower singularity exponent can appear.

Similarly, this study would prot from using very large DNS, giving access to both higher resolutions and higher Reynolds numbers. In order to study such a DNS, it would be necessary to adapt the current post processing code to scale.

Indeed, the resolution of 768 3 used in this thesis is close to the limit for a code parallelized in OpenMP.

Other applications of the Hölder exponent

Turbulence is not the only research eld that can benet from a tool for the detection of singularities. Indeed, such tools based on Hölder regularity have already been used in various domain. In the thesis by Kestener [START_REF] Kestener | Analyse multifractale 2D et 3D à l'aide de la transformation en ondelettes : application en mammographie et en turbulence développée[END_REF], the WTMM method is applied in the context of mammography. Other medical application include the study of human gait [START_REF] Scafetta | Hölder exponent spectra for human gait[END_REF] or arrhythmia [START_REF] Joshi | Arrhythmia classication using local hölder exponents and support vector machine[END_REF].

The method developed during this thesis relies on the multifractal framework and the active volumes, which makes it ne tuned for the case of turbulence.

However, it could be adapted to other adjacent domains, with the example of Magneto Hydro Dynamics (MHD). As we have observed strong events of Hölder exponents at the location of vortex reconnections in Chapter 4, one might expect that a similar tool adapted to MHD would be able to detect events of magnetic reconnection. It is still not known whether solutions to the Navier-Stokes equation can develop singularities from regular initial conditions. In particular, a classical and unsolved problem is to prove that the velocity field is Hölder continuous with some exponent h < 1 (i.e. not necessarily differentiable) at small scales. Different methods have already been proposed to explore the regularity properties of the velocity field, and the estimate of its Hölder exponent h. A first method is to detect of potential singularities via extrema of an "inertial" dissipation D * = lim →0 D I that is independent on viscosity [1]. Another possibility is to use the concept of multifractal analysis that provides fractal dimensions of the subspace of exponents h. However, the multifractal analysis is a global statistical method that only provides a global information about local Hölder exponents, via their probability of occurrence. In order to explore the local regularity properties of a velocity field, we have developed a local statistical analysis, that estimates locally the Hölder continuity. We have compared outcomes of our analysis, with results using the inertial energy dissipation D I . We observe that the dissipation term indeed gets bigger for velocity fields that are less regular according to our estimates. The exact spatial distribution of the local Hölder exponents however shows non trivial behavior and does not exactly match the distribution of the inertial dissipation.

I. INTRODUCTION.

Viscous incompressible fluids are described by the Incompressible Navier-Stokes equations (INSE) in spacetime

∂ t u i + u j ∂ j u i = - 1 ρ ∂ i p + ν∂ j ∂ j u i + f i , (1)
∂ j u j = 0,
where Einstein summation convention over repeated indices is used and u i (x, y, z, t) is the velocity field, p (x, y, z, t) the pressure field, ρ the (constant) mass density, f i (x, y, z, t) some forcing and ν the molecular viscosity. A natural control parameter for the INSE is the Reynolds number Re = LU/ν, which is built using a characteristic length L and velocity U . Mathematically, it is not known whether a solution of the INSE which is smooth at some initial time remains smooth at all later times. This problem was first addressed by Leray [2] who introduced the notion of weak * florian.nguyen@univ-lille.fr † shvydkoy@uic.edu; acheskid@uic.edu solutions (i.e. in the sense of distribution). This notion has since remained a framework to explore regularity condition for both INSE and its inviscid limit, the Euler equation. It prompted people to consider weaker regularity conditions on u. In particular, a classical and unsolved problem is to prove that the velocity field is Hölder continuous with some exponent h < 1 (i.e. not necessarily differentiable) at small scales, i.e. to find under which conditions the following holds:

|u (x + ) -u (x) | < C h . (2)
Let us note that Hölder continuity defined in equation ( 2) is a weaker regularity condition than differentiability and uses the velocity increment δu(x, ) = u(x + )u(x) as a building block. A perhaps more tractable problem is to consider local scaling (Hölder) exponents of the velocity fields. They are defined as:

h (x) = lim →0 ln max | |= |δu (x, ) | ln ( /L) , (3)
where L is a characteristic integral length of scale. Mathematically, these exponents are only defined under the assumption that the limit ν → 0 is taken, before the limit → 0. In the sequel, we consider extension of this definition by intermediate asymptotic in a range of scale

η h(x)
L, where η h(x) ∼ LRe -1/(1+h(x)) . It is however still difficult to estimate directly such local Hölder exponents via a simple fit on velocity increments based on equation (3). Indeed, since accessible resolutions numerically or experimentally do not allow to probe the velocity field at arbitrarily small scale, the range of scales over which the fit is performed is limited. This induces a lot of noise in the determination of h(x), that makes it unreliable.

A global statistical method, was developed by Parisi and Frisch [3] to quantify the probability of observation of a singularity of scaling exponent h via a function C(h). From C(h) we may infer the multifractal spectrum as D(h) = D -C(h), where D is the space dimension. The team of A. Arneodo developed a powerful method based on wavelet transform to measure such spectrum for 1D turbulent velocity fields. This method has been applied to experimental measurements of one velocity component at a single point at high Reynolds numbers in [4], where it was shown that the data are compatible with the multifractal picture, with a most probable h close to 1/3. Later Kestener and Arneodo [5] extended the method to 3D signals (3 components of the velocity field), and showed on a numerical simulation that the picture provided by the 1D measurements was still valid, with the most probable h shifting closer to 1/3. Since the original MFR method is based on statistical properties, it is robust with respect to noise and limited statistics. However, it only provides a global information about local Hölder exponents, via their probability of occurrence. In order to explore the regularity properties of a velocity field, it would however be useful to devise a local statistical analysis, that keeps the robustness of the global MFR method and quantify locally the Hölder continuity, while avoiding pitfalls induced by a naive direct fit of equation ( 3). In the present paper, we extend recent mathematical results of [6] to multifractal settings to provide the best local estimate of Hölder regularity compatible with the global MFR analysis. Our method is in spirit similar to an inference algorithm based on information theory: the local estimate depends on nondimensional constants that are are statistically fitted on the available data set. Therefore, the estimate depends on the data set, and can only provide information a posteriori. The larger the data set, the more information available, the better the estimate. In our case, the larger the data set, the lowest values of h can be explored.

The paper is organized as follows: in Section II, we describe the wavelet based velocity increments and our method for anisotropic multifactal analysis. In Section III, we generalize the notion of active volume introduced by [6] and show how it can be used to build a local estimate of the Hölder exponents . The method is tested using a fractional brownian motion in Section IV. It is further implemented on 3D turbulent velocity fields computed using DNS in Section V. Our results are discussed in Section VI.

II. VELOCITY INCREMENTS AND MULTIFRACTAL

Our method relies heavily on velocity increments and multifractal theory. We therefore provide some generalities about velocity increments in appendix A. We then proceed by defining our new approach using wavelet transforms.

A. Wavelet velocity increments

For an incompressible velocity field, both side of the first equation in (A6) are reduced to zero, so we cannot use it to get any information on φ, the potential part of the velocity field. The usual velocity increments defined in equation (A2) can be seen as the wavelet transform of the velocity with a Haar wavelet. This motivates us to define wavelet velocity increments δW (u) through the smoothed velocity gradient G ij and its symmetric and antisymmetric parts defined as:

G ij = dy 3 ∇ j Ψ y u i (y), S ij = 1 2 G ij + G ji , A ij = 1 2 G ij -G ji , (4)
where Ψ(x) = exp(-x 2 /2)/N is the Gaussian function. By construction, A ij is a function only of the rotational part Q, and S ii is a function only of φ.

We then define normal, longitudinal and transverse wavelet velocity increments as:

δ L W (u) = max ij |S ij |, δ P W (u) = max ij |A ij |, δ G W (u) = δ L W (u) 2 + δ P W (u 2 ) 1/2 . (5)
By construction, δ P W (u) provides information only on the rotational part of the velocity field, while all information regarding the potential part is stored in δ L W (u).

The component δ G W (u) gives general scaling properties about the velocity field. This definition bears some similarity with the definition by Kestener and Arneodo [5], that uses principal values of |G ij |. Like theirs, our definition involves no derivative over velocity fields, so it does not introduce additional noise. However, our definition does not involve computation of singular values, that may generate some noise. So, our definition is the smoothest possible one can imagine that enables to quantify scaling properties of potential and rotational part of a given velocity field.

B. Multifractal theory

The multifractal (hereafter MFR) theory that assumes that the Hölder exponent at scale , defined as

h (x, ) = ln max | |= |δu (x, ) | ln ( /L) , (6)
follows a large deviation property [3] Prob

(h (x, ) = h) ∼ L C(h) . (7)
Formally, C(h) corresponds to the codimension of the set where the local Hölder exponent at scale is equal to h. From C(h) we may infer the MFR spectrum as

D(h) = D -C(h)
, where D is the space dimension.

Computation of the multifractal spectrum

Parisi and Frisch devised in 1987 the first global determination of the multifractal spectrum through the structure functions [3]. The latter are defined from the velocity increments and are supposed to scale as:

|δu| p ∝ ζ(p) . (8)
The set of scaling exponents ζ(p) is then connected to the multifractal spectrum through a Legendre transform:

ζ(p) = min h (ph + D -D(h)), (9)
Assuming the spectrum D(h) is strictly concave, we can apply a variational formula to recover:

h(p) = dζ(p) dp . (10)
Notice that for p = 0, h(0) represents the most probable scaling exponent (with D(h) equals to the full space dimension). When p → -∞, the right part of the MFR spectrum (the more and more regular points) is explored; while when p → +∞, it is the left part of the MFR spectrum (the more and more singular points). This method is global in the sense that it uses global averages of velocity increments to get the multifractal spectrum, and the set of available Hölder exponents. It provides the probability to observe a given Hölder exponent, but provides no information about their exact location. Other methods have been designed based on structure functions. One notable example is the work of Jensen [7] which uses an inverted version of equation ( 8) to obtain potentially more accurate results for p < 0. As we are more interested in probing the less regular region of the velocity field, corresponding to p > 0, we will not use this method in the following.

WTMM method

In 1991, A. Arneodo and his group showed that a much more efficient algorithm of computation of the multifractal spectrum for scalar fields was possible, provided one uses structure functions based on wavelet increments instead of velocity increments. An extension of this method to 3D vector field was developed later by Kestener [5,8] based on singular value decomposition of the wavelet transform of the velocity gradient G ij = ∂ j u i . Here, we apply their algorithm, using the smooth velocity increments equation (5). The algorithm follows the (Wavelet Transform Modulus Maxima) (WTMM) method. We first compute the wavelet transform skeleton, which provides us with a set of lines {L k } where L k ( ) gives the position of the line indexed by k at scale . Then, one computes the partition functions for the different type of wavelet based increments as follow:

Z T (q, ) = L k ∈L( ) δ T W (u)(L k ( )) q (11)
where T stands for G,L or P, q ∈ IR and L( ) is the set of lines that exists at scale .

Those partitions functions behave as a power law of the scale:

Z T (q, ) ∝ τT (q) . (12)
Then by Legendre transforming τ T (q), we get the singularity spectrum D T (h) = min q (qhτ T (q)). Alternatively, we can bypass the Legendre transform by computing:

h T (q, ) = L k ∈L( ) ln |δ T W (u)(L k ( ))|W(q, L k , ), (13)
D T (q, ) = L k ∈L( ) δ T W (u)(L k ( )) ln W(q, L k , ) (14)
where W(q, L k , ) =

1 ZT (q, ) δ T W (u, )(L k ( ))
q is a Boltzmann weight computed from the wavelet transform skeleton. These quantities behave as power laws as follows:

h T (q, ) ∝ hT (q) , (15)
D T (q, ) ∝ DT (q) (16)
and thus we can recover the singularity spectrum D T (h).

III. LOCAL CHARACTERIZATION OF MULTIFRACTAL VECTOR FIELDS

The previous section describes a procedure to compute globally a multifractal spectrum. We now describe a method that gives access to estimates of local scaling properties of the velocity field, using the notion of nested active volume described in appendix B.

A. Active volumes and multifractal

We first define S p,dp as:

S p,dp ( ) = (δW (u)) p+dp (δW (u)) p , 1/dp . ( 17
)
The limit as dp → 0 can be computed as:

S p ( ) ≡ S p,0 ( ) = exp ln δW (u) δW (u) p δW (u) p , . (18)
One can prove that if u is multifractal, then :

ln S p ( ) ≈ h(p) ln , (19)
which gives access to the value h(p).

We can then define the active volume A p by:

x ∈ A p if f δW (u)(x) c p S p ( ), (20)
where c p is a scale independent constant, to be fitted consistently later. The property verified by the points in the volume F h(p) is analogous to the property of equation ( 2). This analogy leads us to define a local singularity exponent h such that all elements x in A p verify h(x) h(p).

As a result, the frontier of this volume, i.e. the isosurface of the increment δW (u) at value c p S p,0 ( ), is identified as the set of the points where the local singularity exponent is h. While the local Hölder exponent h cannot be continuous in the multifractal formalism, the exponent h is continuous by definition. As a result, it is not a real measure of the Hölder exponent. However, it will still be relevant as a measure of the local regularity of the field, hence the name "singularity exponent". This property will be shown by comparing it to other used criteria in the following part.

B. Statistical fitting of the constants

To completely define the field h(x), we need to fix the scale independent constants c p . These constants are the core of the method, as we can tune them such that the field of singularity exponents h provide a meaningful estimation of the real Hölder exponent smoothed at the scale . In order to guarantee that the singularity exponent is physically consistent with the global multifractal properties, we set c p such that the fractal dimensions of the isosurfaces of h match with the multifractal spectrum. More formally:

Dim x| h(x) = h = D(h). (21)
This leads directly to:

Dim ({x|δW (u)(x) = c p S p ( )}) = D(h(p)). (22)
In practice the dimension of the isosurfaces of h, which are also isosurfaces of the velocity increment, can be measured using a box-counting method. If we note D bc (T ) the box-counting dimension of the isosurface of the velocity increment at scale of value T , the problem reduces to the following equation:

D bc (c p S p ( )) = D(h(p)). (23)
The multifractal spectrum are obtained by WTMM, the S p ( ) can be computed from the fields, as well as the function D bc . The quantity h(p) can be accessed using the WTMM or power laws on the S p ( ) as described in equation [START_REF] Duchon | Dissipation d'énergie pour des solutions faibles des équations d'Euler et Navier-Stokes incompressibles[END_REF]. As the two methods are expected to give similar results, we will choose in practice the method of fitting power laws over the functions S p ( ) for simplicity. The functions involved are not monotonous so equation ( 23) could have multiple c p solution. However, we can expect the function D bc (T ) to be concave with a maximum at T max . If we impose c p S p ( ) < T max for p > 0 and T max and c p S p ( ) > T max for p < 0, we get a unique definition of c p . Note that our method ensures that the field h(x) can be recovered without computing explicitly the constant c p . Indeed, we can derive a relationship h(T ) from the following equation:

D bc (T ) = D( h), (24)
which is derived from equation [START_REF] Falconer | Fractal geometry: mathematical foundations and applications[END_REF]. We can ensure that this function is well defined by imposing h(T ) < h(p = 0) for T > T max and reciprocally. While this is the method we will use preferentially to compute maps of h, the knowledge of c p is still pertinent as it is scale independent. This property makes it possible to test the soundness of the method by checking that the procedure, repeated at a different scale, verify the equation [START_REF] Frigo | The design and implementation of FFTW3[END_REF] for the same value of c p . Once a value of c p is computed using one scale, one can then obtain maps of h at any scale in the inertial range for a much lower computational cost. The accuracy of this statements will be discussed in the following.

IV. APPLICATION TO FRACTIONAL BROWNIAN MOTION (FBM)

A. Definition and implementation of fBm

In order to test our analysis, we used synthetic field with well determined fractal behavior. The simplest possible field in this case is the fractional Brownian motion (or fBm), which was already used in [5]. The fBm has properties of isotropy that can be characterized by a single Hölder exponent H for any point.

As a consequence, the multifractal spectrum of a fBm of exponent H in dimension d reduces to a single point:

D(h) = d if h = H, = -∞ otherwise. (25)
In this case, the exponent τ (q) of both the original WTMM method and our variation should have the following expression:

τ (q) = qH -d. (26)
The fBm used in the following are generated by filtering the Fourier transform of a white noise. This method gives a good approximation of the expected invariance properties. For all purposes requiring a range of Hölder exponents, it is assumed that this practical realization with a finite number of snapshots of finite size has a wider multifractal spectrum, centered around the theoretical value H. We chose the value of H = 1/3 for a 2 dimensions, 3 components (2D3C) fBm. The objective in the following is to recover this value H using both the original WTMM and our variation of it.

B. Test of the global MFR analysis

As a first test and benchmark of our methods, we generate a hundred 2D3C fBm (H = 1/3) in squares of size unity with a definition of 4096×4096. The Figure 1 shows the τ (q) corresponding to the different components, as well as the output of the original WTMM method.

The results of our new method are consistent with the theoretical predictions, although a little less accurate than the native method. The consequence of this lack of accuracy is an estimation of the Hölder exponent h ≈ 0.31 using the global term, which is slightly inferior to the value of h ≈ 0.33 evaluated by the original WTMM method. Furthermore, the computation of the perpendicular component suffers from the lack of the third dimension in the provided fields, hence a spectrum spreading over a much wider range of h than expected.

C. Test of the local estimate

In parallel, we can apply our second method to extract the local coefficients h in one field of the same set of fBm. The computation requires to chose a value of . In the sequel, we fix the value = 2.8 × 10 -3 for any computation on the fBm at definite . This value, corresponding to 11.4 grid points at the resolution of the field, is small enough to probe the small structures and large enough so that wavelets are well resolved. The Figure 2 shows the function D bc (T ) resulting from the boxcounting, while the Figure 3 illustrates the power law behavior of S p ( ). The values of h are convincingly independent with respect to p, with very little fluctuations. The value of h obtained is around h ≈ 0.34, which is slightly higher than the expected value of H = 1/3 but still within acceptable range.

Applying the method described to compute fields of h would normally require to match the box-counting dimensions from Figure 2 with the multifractal spectrum Exponent τ as a function of the parameter q for the global (red circles), longitudinal (green squares), perpendicular (blue triangles) components as well as for the native method (black diamonds) applied to 100 fields of fBm in 2 dimensions with 3 components of velocity at H = 1/3. The 3 components are generated independently before enforcing the divergence free condition. The theoretical curve from equation ( 26) is materialized by the solid magenta line. Inset figure: Corresponding parametric plot of D(q) with respect to h(q). The dashed line materializes the theoretical value of H = 1/3. from Figure 1. Unfortunately, the degenerated spectrum makes it impractical. A solution is to use the fact that we are working with a monofractal and ensure that for all p, c p S p ( ) = T max . Considering the fact that it is very difficult to get measurements of box counting dimension that reach 2, the dimension of the space, we also take the liberty to renormalize D bc (T ) such that D bc (T max ) = 2.

In this artificial case, analyzing maps of h will not be very instructive, so we verify the validity of using the scale invariant coefficients c p in order to compute fields of h at a different scale. To such aim, we use instead the box-counting dimension at the lower scale 2 = 8.2×10 -4 corresponding to approximately 3.4 grid points. This allows to check that equation ( 23) is satisfied, thus proving the point in the case of fractional Brownian motion.

The Figure 4 shows the box-counting dimension of the isosurfaces of "velocity" increments computed at the scale 2 . Although the quality of the measure is lower for less resolved wavelets, we can still verify that the pseudo multifractal spectra constructed using the scaling exponent from the S p ( ), the c p obtained at higher scale and the box counting dimensions from Figure 4, is coherent with the theoretical spectrum, as illustrated in Figure 5. Now that the methods have been tested in the case of fBms, we apply them in the more complex and realistic case of turbulence fields. The red circles (resp. green squares, blue triangles) correspond to the spectrum reconstructed using the global (resp. longitudinal, perpendicular) increments. The solid vertical line materializes the theoretical h = 1/3.

V. APPLICATION TO VELOCITY FIELDS IN NUMERICAL SIMULATION

A. Characteristics of the numerical simulation

We generate turbulent velocity fields via two simulations of turbulence in a box of size 2π. The corresponding parameters are presented in Table I. We use a pseudo spectral code [9] with a forcing term with the symmetries of the Taylor-Green vortex. The first simulation (Run I) provides velocity fields with an inertial range large enough to validate our method against known scaling laws. The second simulation (Run II) is very well resolved with an extensive dissipative range. It will allow us to probe the actual singular behavior of the flow at small scales of the order of the Kolmogorov scale. The energy spectra presented Figure 6 illustrate the difference of scale resolution between these two simulations.

B. MFR spectra

By applying the previously described variation of the WTMM method, we extract the global multifractal singularity spectra out of 10 snapshots of velocity fields for the two simulations, presented in Figure 7 and8. As the velocity fields include 768 3 grid points and are decorrelated in time, this represents enough statistics to guarantee the convergence of the following results.

In the inertial regime, we recover spectra centered around h = 0.40 for all components but the perpendicular one, which is higher than the expected value of h = 1/3 according to the Kolmogorov 41 theory. The analysis by Arneodo [10] The vertical line materializes the expected exponent of h = 1/3 according to K41. The error bars correspond to a shift of the fitting range by 5% for the power laws. Inset: Exponents τ as a function of q for the global (red circles), longitudinal (green squares)and perpendicular (blue triangle) increments, as well as the native method (black diamonds). The error bars correspond to a shift of the fitting range by 10% for the power laws. Inset: Exponents τ as a function of q for the global (red circles), longitudinal (green squares)and perpendicular (blue triangle) increments, as well as the native method (black diamonds).

exponent of around h = 0.34, which was much closer to the Kolmogorov value. It is however not clear how universal this value should be, in particular as we use data from anisotropic turbulence where Arneodo used homogeneous turbulence. The fact that the perpendicular component does not follow exactly the same tendency is the manifestation of the anisotropy of the flow at this scale simulated with a Taylor-Green forcing. Additionally, one can also notice the remarkable agreement between the native method from Pierre Kestener and the modified method applied to global increments, which further validates the choice of the increments in equation (5).

The multifractal spectra for the simulation resolved under the Kolmogorov scale provide significantly different results, shown Figure 8. Indeed, the values of the singularity exponents are much higher since the velocity field is much more regular at this scale. As a reminder, a Hölder exponent 2 < h < 3 means that the velocity field at this point is twice differentiable but not thrice. We also observe that the spectra reach a maximum D ≈ 1.4 < 3. This can be interpreted as singularities appearing less likely once the behavior of the velocity fields at the dissipative scale is taken into account. Finally, the superposition of the spectra for the different components hint the recovery of isotropy at small scale.

C. Estimates of local singularity exponents

As already discussed, the multifractal spectra obtained here do not provide any local information. While the field h defined with equation ( 21) is not equivalent to a measure of local Hölder exponents, we will show in the following that it is still relevant for the measure of local regularity.

To compute this field, we follow the same steps as for the fBms, and first compute boxcounting dimensions over isosurfaces of the increments (Figure 9), as well as the power law behavior of the coefficients S p ( ) (Figure 10). The computation of increments before applying the boxcounting algorithm implies to chose a value of . The value chosen in the following is ≈ 0.24, well within the inertial range.

The situation here slightly differs from the fBm case, since the multifractal spectrum is non degenerated. This allows us to match the dimensions from figures 9 and 7 using equation [START_REF] Frisch | Singularities of Euler ow? Not out of the blue![END_REF] to get a relationship directly connecting the velocity increments to the coefficient h. By taking care of associating the high velocity increments to low values of h and reciprocally, we obtain the relationship represented Figure 11. At this point, we do not produce separate measurements for the perpendicular and longitudinal components anymore, as they only provide partial information about the velocity field.

Using this procedure, we are able to deduce maps of h from the corresponding maps of velocity increments. One flaw of the process however is that the range of T shown in Figure 11 is smaller than the total range of velocity increments encountered in practice. The reason for this comes from the fact that both WTMM and the boxcounting of isosurfaces of the velocity increments are limited by the statistics available. The lower the fractal dimension, the harder it is to get enough statistics. With the amount of statistics available (10 fields of size 768 3 ), we cannot get reliable value of boxcounting dimension below approximately 2.4 as shown in Figure 9. As a result, we cannot get values of h for velocity increments much higher than 0.5. This will eventually leads to the maps of singularity exponents having regions where we know that h < 0.1 but cannot give any precise value.

D. About the scale invariance of the cp coefficients

We have stated in Section III that the coefficients c p defined in equation ( 20 12. Pseudo spectra reconstituted using the coefficient cp computed at scale ≈ 0.24 and the boxcounting dimensions of the velocity increments computed at scales ≈ 0.12 (dashed green curve), ≈ 0.24 (solid red curve) and ≈ 0.46 (dotted blue curve). The scale are also materialized by vertical lines of the same color on the Figure 10.

sured at a given scale to extrapolate to all scales, thus saving computation time. Indeed, once these coefficients are known, the knowledge of the S p ( ) is sufficient to recover the active volume at any scale without having to compute box-counting dimensions nor applying the WTMM method.

We use the present simulation data to test the validity of this property. The first hurdle here is that the coefficient h(p) from equation [START_REF] Duchon | Dissipation d'énergie pour des solutions faibles des équations d'Euler et Navier-Stokes incompressibles[END_REF] as well as the multifractal spectrum from the WTMM method are defined in sections II and III as dependent of the scale . However, both actually requires the estimation of the power law behavior with respects to , which are done in practice by fitting curves over a range of . The limited accuracy over the power law exponents, as well as the actual corresponding scale within the fitting range used, are expected to reduce the accuracy of the c p computed.

As the c p do not have an obvious physical meaning, it is difficult to estimate what range of error in the measures of c p would be coherent with the property of scale independence. This leads us to use a slightly different method to test this property. Using the results computed at the scale of ≈ 0.24 used in the previous figures, we compute the coefficients c p that we assume to be scale independent. Using these coefficients and the box-counting dimensions of the isosurfaces of velocity increments computed at two different scales, we recover pseudo multifractal spectra using equation [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. The results are presented in Figure 12. We observe that the spectra are qualitatively coherent. In particular, the green curve shows that when we reach the lower bound of the inertial range, the spectrum widen with much higher values of h. In other words, the velocity fields starts to get regularized by viscosity.

A direct comparison with the output of the WTMM method is unfortunately not possible, as the WTMM performs poorly with fitting ranges too short. More statistics would be required to perform the power law fit of equations ( 15) and ( 16) with a good accuracy on a very small fitting range. The result of this analysis is that using coefficient c p computed at a different scale will provide qualitatively sound results. Since the complete method is dependent on fitting scaling exponents, it is reasonable to use the c p for quantitative measurements.

E. Comparison with another indicator of regularity

By construction, our method guarantees that the local estimates of the Hölder exponent are statistically meaningful and respect the global MFR properties. To evaluate the physical soundness of our estimate (i.e. how good it is to detect area of lesser regularity), we may compare our maps of h(x) with maps of local energy transfer [1], hereafter referred to as Duchon-Robert term. Indeed, this quantity has been used in [11] to build a new criterion to detect areas with Hölder exponent h < 1/3, by looking at local maxima of such quantities. The Duchon-Robert term is expressed as:

D I (x) = 1 4 dy∇Ψ (y) • δu(x, y)||δu(x, y)|| 2 (27)
where Ψ (x) = Ψ(x/ )/ 3 and δu(x, y) is the conventional velocity increment as defined in equation A2.

When goes to 0, the term D I (x) scales as 3h-1 , hence it diverges for an Hölder exponent h < 1/3. In practice, at any given scale , the higher the term D I (x) the less regular the velocity field is in x. So if our new index h is meaningful, its areas of lower values should correlated with areas of higher value for the Duchon-Robert term.

We computed maps of both the singularity exponent h and the Duchon-Robert term D I (x). The Figure 13 illustrates one of the strong events that can be found in our velocity fields. It appears in this case that the two quantities are correlated, meaning that regions of lower h visually corresponds to regions of high |D I (x)| in absolute value. However, looking in more detail, the region of low h seems to be enclosed between two regions of high |D I (x)|.

While it appears that the minima of h are not located at the same location as the extrema of |D I (x)|, it is still relevant to observe the joint PDF of those two terms in order to confirm this observation. The results of statistics aggregated over 50 fields extracted at regular interval over approximatively 30 turnover times are reported in Figure 14. We estimate that the small scales of two different fields are not correlated. While the conditional PDF is not perfectly converged for the extreme values of h, it is enough to do qualitative observations. We observe a very clear correlation between the two quantities. In particular, lower values of h correspond to higher values of |D I (x)| in absolute value, which can be associated to a more singular velocity field. This observation comforts us in the idea that the field of singularity exponents h can be used as a measure of the local regularity.

VI. DISCUSSION

In this paper, we have introduced a tool to estimate the regularity of velocity fields. It is derived using the Wavelet Transform Modulus Maxima method and associates the multifractal spectra and the isosurfaces of velocity increments to compute a field h(x) which has similar properties than the Hölder exponents. We have checked the physical soundness of our estimate by comparing it with the energy transfer term |D I (x)| [1], that can be used as a marker of regions of lesser regularity [11]. We have found that our h(x) and |D I (x)| are globally statistically correlated, but that local maxima of |D I (x)| do not coincide exactly with local minima of h(x). On the practical side, we found that, once the relation between the wavelet based velocity increments and the coefficients h shown Figure 11 is established, computing h(x) from a new velocity field is approximately three times faster than computing the field of D I (x). As a result, our method can be used to process velocity fields for candidates of singularities faster than through the energy transfer term.

The study of particular events like the one presented in Figure 13 shows a non trivial correlation of the low h regions with respect to the extreme D I regions. This encourages further studies using both scalars to characterize the velocity structures and the regularity properties of turbulent flows. Finally, this study has been focused on turbulence driven by a Taylor-Green forcing. It could be interesting to study other types of flows like turbulence induced by Rayleigh-Taylor and Richtmyer-Meshkov in-stabilities [12,13]. Using this and that on Ω T \A , |f | α we obtain

Lemma B.1. There exists a set

A ⊂ Ω T with |A| = V r,p |Ω T | such that (1 -c r,p ) ΩT |f | p dµ A |f | p dµ, (
|f | p χ ΩT \A α p-r |f | r χ ΩT \A 1 V p-r r r,p |f | r χ A p-r r |f | r χ ΩT \A |f | p |f | r χ ΩT \A r p |f | r χ A p-r r |f | r
Note that the latter fraction is of the form θ r p (1θ) p-r r , θ ∈ [0, 1], which attains its maximum exactly at the value c r,p . This proves the lemma.

Note that as r → p, c r,p → 1, and as a result the information about concentration of f gets lost. However, the corresponding concentration factors converge to something non-trivial, namely,

V p = lim r→p V r,p = |f | p exp - |f | p ln |f | p |f | p . ( B4
)
Let us call it p-concentration factor. The factor have some natural monotonicity properties, which will be addressed in the next section.

Amplitude

Although Lemma B.1 explains why the measure of concentration of f can be defined by V r,p , the method of proof provides little constructive information about the threshold amplitude α. Such amplitude can be defined quite explicitly, if one is ready to sacrifice precise measure of the set A. Namely, let us denote

s r,p = |f | p 1 p-r |f | r 1 p-r . ( B5
)
We will call it (r, p)-amplitude. Note that the physical units of s r,p and f coincide.

Let us further expand our dictionary by calling a point (x, t) ∈ Ω T (r, p)-active if the amplitude of f passes the s r,p threshold: |f (x, t)| > c r,p s r,p , where 0 < c r,p < 1 are empirical adimensional factors. Collectively the set of all (r, p)-active points form a (r, p)active domain: (B9)

The amplitudes have several monotonicity properties. First, applying Jensen's inequality,

s p-r r = exp |f | r ln |f | p-r |f | r |f | r exp(ln |f | p-r ) |f | r = s p-r r,p .
Thus, s r s r,p , for all r < p. Second, one versifies again by Jensen, that ∂ r s r,p , ∂ p s r,p 0. By taking the limit as r → p, s r,p s p . This in turn implies that s r s r,p s p , for all r < p.

From (B9) we see that ln s r,p is precisely the mean value of ln s τ on [r, p]. In view of the identities

s p r,p V r,p = |f | p , s r r,p V r,p = |f | r
we deduce the opposite monotonicity properties of the concentration factors: ∂ r V r,p , ∂ p V r,p 0, and thus V r V r,p V p , for all r < p.

p-active regions

Let U 0 be a characteristic size of f , and let us assume that the units of f is velocity as well. It is then more natural to rewrite the formula for s p as follows:

s p = U 0 exp    |f | p ln |f | U0 |f | p    . (B10)
(note that this holds for any U 0 > 0). This way of expressing s p makes the main exponent adimensional. We can view the ratio s p /U 0 as an adimensional threshold value for an active size. Next, we show the active region defined in [START_REF] Grafke | Numerical simulations of possible nite time singularities in the incompressible Euler equations: comparison of numerical methods[END_REF]

I. INTRODUCTION

A well-known feature of any turbulent flow is the Kolmogorov-Richardson cascade by which energy is transferred from large to small length scales until the Kolmogorov length scales below which it is removed by viscous dissipation. This energy cascade is a non-linear and an out-of-equilibrium universal process. Moreover, the corresponding non-dimensional energy spectrum E(k)/ǫ 2/3 η 5/3

is an universal function of kη, where η = (ν 3 /ǫ) 1/4 is the Kolmogorov length scale, ǫ the mean energy dissipation rate per unit mass, and ν the kinematic viscosity. However, there seem to be little dependences on the Reynolds number, boundary, isotropy or homogeneity conditions 2 . In facts, the energy spectrum is based upon a quantity, the velocity correlation, that is quadratic in velocity. Nevertheless, it is now well admitted that the universality does not carry over for statistical quantities that involve higher order moments. For example, the velocity structure functions of order p, given by S p (ℓ) = u(x + r)-u(x) p x, r =ℓ are not universal, at least when expressed in units of the Komogorov scale η and velocity u K = (νǫ) 1/4 (see below, section III B for an illustration).

The mechanism behind this universality breaking was identified by, wherein a generalization of the Kolmogorov theory introduced, based on the hypothesis that a turbulent flow is multifractal.

In this model, the velocity field is characterized locally by an exponent h, such that |δ ℓ u(x)| ≡ u(x + r)u(x) r =ℓ ∼ ℓ h(x) ; here h is a stochastic function that follows a large deviation property 3 P (log(|δ ℓ u|/u 0 ) = h log (ℓ/L 0 )) ∼ (ℓ/L 0 ) C(h) , where u 0 (resp. L 0 ) is the caracteristic integral velocity (resp. length), and C(h) is the multifractal spectrum. Velocity fields with h < 1 are rough in the limit ℓ → 0. In real flows, any rough field with h > -1 can be regularized at sufficiently small scale (the "viscous scale") by viscosity. The first computation of such dissipative scale was performed by Paladin and Vulpiani, who showed that it scales with viscosity like η h ∼ ν 1/(1+h) , thereby generalizing the Kolmogorov scale, which corresponds to h = 1/3. Such a dissipative scale fluctuates in space and time (along with h), resulting in non-universality for high order moments, at least when expressed in units of η and u K . further suggested that the prediction can be extended to the velocity structure functions, so that, at any given p log(S p (ℓ)/u p K )/ log(Re) should be a universal function of log(ℓ/η)/ log(Re); they found good agreement for p up to 6. The velocity measurements, in the above experiments, were performed using hot wire anemometry, which provide access to only one component of velocity.

To our knowledge, no further attempts have been made to check the claim with more realistic measurements.

The purpose of the present paper is to reexamine this claim; however, now using the velocity fields obtained from the Stereoscopic Particule Image Velocimetry (SPIV) in experiments and the direct numerical simulations (DNS) of the Navier-Stokes equations (NSE). We show that the numerical and experimental velocity structure functions up to order 9 follow a log-universality 1 ; they indeed collapse on a universal curve, if we use units that include log(Re) dependence. We then investigate the meaning and consequences of such a log-universality, and show that it is connected with the properties of a "multifractal free energy", based on an analogy between multifractal and thermodynamics (see for summary). We show that in such a framework, the existence of a fluctuating dissipation length scale is associated with a phase transition describing the relaminarisation of rough velocity fields with different Hölder exponents.

II. EXPERIMENTAL AND NUMERICAL SETUP A. Experimental facilities and parameters

We use experimental velocity field described in 5 . The radial, axial and azimuthal velocity are measured in a von Krmn flow, using Stereoscopic Particule Image Velocimetry technique at different resolutions ∆x. The von Krmn flow is generated in a cylindrical tank of radius R = 10 cm through counter-rotation of two independent impellers with curved blades. The flow was maintained in a turbulent state at high Reynolds number by two independent impellers, rotating at frequency F . Figure 1 shows the sketch of the experimental setup. The five experiments were performed in conditions so that the non-dimensional mean energy dissipation per unit mass is constant. The viscosity was monitored using mixture of water and glycerol, so as to vary the Kolmogorov length η. ∇u 2 , the mean squared velocity u rms and the viscosity ν. All velocity measurements are performed in a vertical plane that contains the rotation axis. The case A corresponds to measurements over the whole plane contained in between the two impellers, and extending from one side to the other side of the cylinder. Its resolution is 5 to 10 times coarser than similar measurements performed by zooming on a region centered around the symmetry point of the experiment (on the rotation axis, half way in between the two impellers), over a square window of size 2 cm × 2 cm.

B. Direct Numerical Simulation

The direct numerical simulations (DNS), based on pseudo-spectral method, were performed in order to compare with our experimental data. The DNS runs with R λ = 25, R λ = 80, R λ = and R λ = 138 were performed using the NSE solver VIKSHOBHA was carried out using another independent pseudo-spectral method based NSE solver. The velocity field u was computed on a 2π triply-periodic box.

Turbulent flow in a statistically steady state was obtained by using the Taylor-Green type external forcing in the NSE at wavenumber k f = 1 and amplitude f 0 = 0.12, the value of viscosity was varied in order to obtain different values of R λ (see Ref. 6 for more details).

III. THEORETICAL BACKGROUND

A. Velocity increments vs Wavelet Transform (WT) of velocity gradients

The classical theories of Kolmogorov 7? are based on the scaling properties of the velocity increment, defined as δ ℓ u = u(x + r)u(x), where ℓ = |r| is the distance over which the increment is taken.

As pointed out by, a more natural tool to characterize the local scaling properties of the velocity field is the wavelet transform of the tensor ∂ j u i , defined as:

G ij (x, ℓ) = dr∇ j Φ ℓ (r) u i (x + r), (1)
where Φ ℓ (x) = ℓ -3 Φ(x/ℓ) is a smooth function, non-negative with unit integral. In what follows, we choose a Gaussian function Φ(x) = exp(-x 2 /2)/(2π) 3 2 such that Φ(r)dr = 1. We then compute the wavelet velocity increments as

δW (u)(x, ℓ) = ℓ max ij |G ij (x, ℓ)|. ( 2
)
This formulation is especially well suited for the analysis of the experimental velocity field, as it naturally allows to average out the noise. It has been verified that the wavelet based approach yields the same values for the scaling exponents as those computed from the velocity increments 6 .

B. K41 and K62 universality

In the first theory of Kolmogorov 7 , the turbulence properties depend only on two parameters: the mean energy dissipation per unit mass ǫ and the viscosity ν . The only velocity and length unit that one can build using these quantities are the Kolmogorov length η = (ν 3 /ǫ) 1/4 and velocity u K = (ǫν) 1/4 . The structure functions are then self-similar in the inertial range η ≪ ℓ ≪ L 0 , where L 0 is the integral scale, and follow the universal scalings:

S p (ℓ) ≡ |δ ℓ u| p ∼ u p K ℓ η p/3 , (3)
which can also be recast into:

Sp ≡ S p S p/3 3 = C p , (4)
where C p is a (non universal) constant.

where τ (p) = µp(3p) and A p is a constant. Such a formulation already predicts an interesting universality, if L = L 0 , as we should have:

L 0 η τ (p) Sp ∼ A p ℓ η τ (p) . (7)
Therefore, we should be able to collapse all structure functions, at different Reynolds number by plotting ( L 0 η ) τ (p) Sp as a function of ℓ η , given that L 0 /η ∼ Re 3/4 . There is however no clear prediction about the value of L and we show in the data analysis section that L differs from L 0 .

The relation (7) shows that log L 0 η τ (p) Sp is a linear function of log( ℓ η ). In principle, such universal scaling is not valid outside the inertial range, i.e. for example when ℓ < η. To be more general than previously thought, it can however be shown using the multifractal formalism as first shown by 4 .

C. Multifractal and fluctuating dissipation length

For the multifractal (MFR) model, it is assumed that the turbulence is locally self-similar, so that there exists a scalar field h(x, ℓ, t), such that h (x, t, ℓ) = log (|δ ℓ u(x, t)|/u 0 ) log(ℓ/L) ,

for a range of scales in a suitable "inertial range" η h ≪ ℓ ≪ L, where L is a characteristic integral-length-scale, η h a cut-off length scale, and u 0 a characteristic large-scale velocity. This scale is a generalization of the Kolmogorov scale, and is defined as the scale where the local Reynolds number ℓ|δ ℓ u|/ν is equal to 1. Writing δ ℓ u = |δ ℓ u| = u 0 (ℓ/L) h leads to the expression of η h as a function of the global Reynolds number Re = u 0 L/ν as η h ∼ LRe -1/(1+h) . This scale thus appears as a fluctuating cut-off which depends on the scaling exponent and therefore on x.

This is the generalization of the Kolmogorov scale η ∼ ν 3/4 ≡ η 1 3

, and was first proposed in.

Below η h , the velocity field becomes laminar, and |δ ℓ u| ∝ ℓ. When the velocity field is turbulent, h ≡ log(|δ ℓ u|/u 0 )/ log(ℓ/L) varies stochastically as a function of space and time. Also, if the turbulence is statistically homogeneous, stationary and isotropic, h only depends on ℓ, the scale magnitude. Therefore, formally, h can be regarded as a continuous stochastic process labeled by log(ℓ/L). By Kramer's theorem 8 , one sees that as in the limit ℓ → 0, log(L/ℓ) → ∞, we have

P [log(δ ℓ u/u 0 ) = h log(ℓ/L)] ∼ e log(ℓ/L)C(h) = ℓ L C(h) , (9)
where C(h) is the rate function of h, also called multifractal spectrum. Formally, C(h) can be interpreted as the co-dimension of the set where the local Hölder exponent at scale ℓ is equal to h.

Using Gärtner-Elis theorem 8 , one can connect C and the velocity structure functions as:

S p (ℓ) = (δ ℓ u) p = hmax h min u p 0 ℓ L ph+C(h) dh. ( 10
)
To proceed further and make connection with previous section, we set ǫ = u 3 0 /L so that S p (ℓ) can now be written:

To understand the nature of the correction, we can compute the value of h such that ℓ = η(h). It

where τ (p, θ) = τ p when θ ≤ 1/(1+h max ) and τ (p, θ) = p(θ-1/3)+C(-1+1/θ) when 1/(1+h max ) ≤ θ ≤ 1/(1 + h min ). As discussed by 4 , this implies a new form of universality that extends beyond the inertial range, into the so-called extended dissipative range, as; log( Sp ) log(Re) = -τ (p, θ)θ, θ = log(L/ℓ)/ log(Re). ( 14)

If the scale L is constant and equal to L 0 , the integral scale, then we have Re = (L 0 /η) 4/3 and the multifractal universality implies that log( Sp )/ log(L 0 /η) is a function of log(ℓ/η)/ log(L 0 /η).

When the function is linear, we thus recover the K62 universality. The multifractal universaility is thus a generalization of the K62 universality.

This form of universality is however not easy to test, as the scale L is not known a priori, and may still depend on Re. In what follows, we demonstrate a new form of universality, that allows more freedom upon L and encompass both K62 and multifractal universality.

D. General universality

Using the hypothesis that turbulence maximizes some energy transfer in the scale space, Castaing 9

suggested a new form of universality for the structure functions, that reads:

where S 0p are some constants and β and F are general functions. Comparing with the K62 or MFR universality formulae (7) or (14), we see that formula ( 16) is a generalization of these two universality with L = L 0 . It allows however more flexibility than K62 or MFR universality through the function β(Re), that is a new fitting function. We test these predictions in Section IV and provide a physical interpretation of ( 16) in Section V.

IV. CHECK OF UNIVERSALITY USING DATA ANALYSIS

The various universality are tested using the velocity structure functions based on the wavelet velocity increments Eq. ( 2), in order to minimize the noise in the experimental data. We define:

S p = |δW (u)(x, ℓ)| p . ( 17
)
We then apply this formula to both experimental data (Table I) and numerical data (Table II), to get wavelet velocity structure functions at various scales and Reynolds numbers.

This is shown in figure 2 for both experimental and numerical data. Obviously, the data do not collapse on a universal curve, meaning that K41 universality does not hold. This is well known, and is connected to intermittency effects 10 .

B. Check of K62 universality

The K62 universality (7) can be checked by plotting:

log L 0 η τ (p) Sp = F log ℓ η . (19)
The collapse depends directly on τ (p), the intermittency exponents. Obtaining the best collapse of all curves is in fact a way to fit the best scaling exponents τ (p). We thus implemented a minimization algorithm that provides the values of τ (p) that minimized the distance between the III. The best collapse is shown on Figure 3-a for the DNS, and Figure 3-b for the experiment. The collapse is better for experiments than for the DNS. However, in both cases, there are significant differences in between points at different R λ , at larger scales, showing that universality is not yet reached.

C. Check of General Universality

We can now check the most general universality, by plotting:

β(Re) log( Sp /S 0p ) log(L 0 /η) = F p, β(Re) log(ℓ/η) log(L 0 /η) , (20)
In this case, best collapse is obtained by fitting two families of parameters: S 0p , β(Re) that were obtained through a procedure of minimization. We take the DNS at R λ = 138 as the reference case, and find for both DNS and experiments, the values of β(Re) and S 0p that best collapse the curves.

The corresponding collapses are provided in figure 5. The collapse is good for any value of Re, except for the DNS at the lowest Reynolds number, which does not collapse in the far dissipative range. The curve look like the portion of a parabola, corresponding to a log-normal statistics for the wavelet velocity increments. Specifically, fitting by the shape:

C(h) = (h -a) 2 2b , (23)
we get a = 0.35 and b = 0.045. This parabolic fit also provides a good fit of the scaling exponents, as shown in figure 4 by performing Legendre transform of C(h) given by equation [START_REF] Frisch | Turbulence: the legacy of AN Kolmogorov[END_REF].

V. THERMODYNAMICS AND TURBULENCE

A. Thermodynamical analogy Multifractals obey a well-known thermodynamical analogy 12? ,13 that will be useful to interpret and extend the general universality unraveled in the previous section. Indeed, consider the quantity:

µ ℓ = |δW ℓ | 3 |δW ℓ | 3 . ( 24
)
By definition µ ℓ is positive definite and µ ℓ = 1 for any ℓ. It is therefore the meaning of a scale dependent measure. It then also follows a large-deviation property as:

P [log(µ ℓ ) = E log(ℓ/η)] ∼ e log(ℓ/η)S(E) ,

so that Z = S3p is the partition function associated to the variable log(µ ℓ ), at the pseudo-inverse temperature p = 1/kT . Taking the logarithm of Z, we then get the free energy F as:

log( S3p ) = F. (27)
By the Gärtner-Elis theorem, F is the Legendre transform of the energy S: F = min E (pE -S(E)).

The free energy a priori depends on the temperature; i.e. on T = 1/kp, on the volume V = log(ℓ/η) and on the number of degrees of freedom system N . If we identify N = (1/β(Re)) log(L 0 /η), we see that the general universality means:

F (T, V, N ) = N F (T, V N , 1), (28)
i.e. can be interpreted as extensivity of the free energy.

Figure 1 . 1 :

 11 Figure 1.1: Schematic illustration of the energy spectra in a turbulent ow in log-log.

Figure 1 . 2 :

 12 Figure 1.2: Schematic of anomalous dissipation scaling.

Figure 2 . 1 :

 21 Figure 2.1: Schematic representation of a VonKármán setup. Figure reproduced from [40].

Figure 2 . 2 :

 22 Figure 2.2: Four types of singularities found among the extreme events extracted from the data at the dissipative scale, where the Kolmogorov scale is resolved. Main gures: colormap of the dimensionless D I (u) eld. The arrows correspond to in plane velocity. The white square materializes the location with the largest turbulent force in norm. Vector eld insets: Velocity eld near the singularity. The arrows stand for the radial (x) and axial (z) components while the colour scale stand for the azimuthal (y) component. Plot insets: D I (u) in dimensionless unit near the singularity as a function of /η.Figure reproduced from [41].

Figure 2 . 3 :

 23 Figure 2.3: Four types of singularities found among the extreme events extracted from the data at the dissipative scale, where the Kolmogorov scale is resolved. Main gures: dimensionless F vector eld for ≈ 5η. The arrows stand for the radial (x) and axial (z) components while the colour scale stand for the azimuthal (y) component. The white square materialize the location with the largest turbulent force in norm. Vector eld insets: Velocity eld near the singularity with the same conventions. Plot insets: Maximum of F in dimensionless unit near the singularity as a function of /η.

Figure 2 . 4 :

 24 Figure 2.4: QR plane with the dierent topologies and three Vieillefosse lines [71] drawn.

Figure 3 . 1 :

 31 Figure 3.1: Two functions with h = 0.5 at x = 0.

Figure 3 . 2 :

 32 Figure 3.2: Representation of the two wavelets in 1D.

Figure 3 . 3 :

 33 Figure 3.3: Three steps of box counting applied to the Sierpinski triangle (drawn here with 5 iterations). If the main triangle is of size 30, the box sizes are from left to right 5, 2.4 and 1.2. The respecting box counts are 24, 76 and 240.

(3. 43 )

 43 where D bc (T ) is the box-counting dimension of the isosurface of the velocity increment at scale of value δW (u) = T .

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: Main gure: Exponent τ as a function of the parameter q for the global (red), longitudinal (green), perpendicular (blue) components as well as for the original WTMM method (black) applied to 100 elds of fBm in 2 dimensions with 3 components of velocity at H = 1/3. The 3 components are generated independently before enforcing the divergence free condition. The theoretical curve from equation (3.56) is materialized in magenta. Inset gure: Corresponding parametric plot of D(q) with respect to h(q). The dashed line materializes the theoretical value of H = 1/3.

Figure 3 . 6 :

 36 Figure 3.6: Values of the threshold T p ( ) for the global increments as a function of scale for several p for the fBm.

Figure 3 . 7 :Figure 3 . 8 :

 3738 Figure 3.7: Box counting dimension of the isocontours of the dierent components of the wavelet based velocity increments for = 8.2 × 10 -4 for the fBm. The color code for the dierent components is the same as for gure 3.6.

Figure 4 . 1 :

 41 Figure 4.1: Time evolution of the energy dissipation rate as a function of time for the Run I. The vertical line materializes the time at which the rst snapshot of velocity was extracted.

Figure 4 . 2 :

 42 Figure 4.2: Energy spectra for the two simulations. The blue continuous (resp. red dashed) curve stands for the simulation in the inertial (resp. dissipative) range. The black line materializes the k -5/3 slope for the inertial range. The horizontal lines of the same color as the curves materialize the corresponding tting range used when computing power laws. The vertical lines materialize the Taylor scale wave number k λ for each simulation.

Figure 4 . 3 :

 43 Figure 4.3: Main gure: Multifractal spectra for the velocity elds from simulation in the inertial region. The red (resp. green, blue) curve corresponds to the global (resp. longitudinal, perpendicular) increments. The black curve corresponds to the increments based on singular value decomposition used in the original WTMM method. The vertical line materializes the expected exponent of h = 1/3 according to K41. The error bars correspond to a shift of the tting range by 5% for the power laws. Inset: Exponents τ as a function of q.

3 .

 3 The power law tting required from equations (3.26) and (3.27) (as well as equation (3.22) ) are performed over the range of scale specied in gure 4.2, well within the inertial range. While a Gaussian wavelet might work in this case as we expect h < 1, we use a Mexican wavelet for the sake of comparison with the dissipative range.

Figure 4 . 4 :

 44 Figure 4.4: (a) Box count of 6 isosurfaces of velocity increments as a function of the box size. The results are totaled over 10 snapshots of velocity. (b) Box counting dimension of the isosurfaces of the dierent components of the wavelet based velocity increments for ≈ 0.24 for the simulation in the inertial range. The color code for the dierent components is the same as for the previous gures.

Figure 4 Figure 4 . 6 :

 446 Figure 4.5: (a) Values of the threshold T p ( ) for the global increments as a function of scale for several p computed over 10 snapshots of run I. The thresholds are renormalized by T 0 ( ). Fits corresponding to the scales delimited by the green and blue vertical lines are given in dotted lines. (b) Reference thresholds T 0 ( ). The t used in the calibration of the method is materialized by a dotted line.

Figure 4 . 7 :

 47 Figure 4.7: Region with strong irregularity extracted from one snapshot of Run I. (a) Velocity along a slice passing through the point of maximal velocity increment (corresponding to an undetermined minimal value of Hölder exponent). The arrows stand for the velocity in plane and the colormap for the third component of velocity. (b) Local energy transfers on the same slice. (c) Local Hölder exponents on the same slice. White area corresponds to velocity increments so strong we cannot associate a value of h using the mapping function of gure 4.6. (d) 3D representation of the irregular region. The magenta isosurface corresponds to h = 0.25. The red isosurface corresponds to D I = 4 .

Figure 4 . 8 :

 48 Figure 4.8: Joint PDF of h and the DuchonRobert energy transfer term D I over 50 elds regularly spaced over approximately 30 turnover times. In order to enlighten the correlation of those terms independently of the number of occurrence, the PDF is rescaled such that every vertical line is a conditional PDF of D I at given h.

(4. 3 )

 3 with Ψ mex the 3D Mexican wavelet at scale . The choice of a Mexican wavelet instead of a Gaussian wavelet is made in preparation for the study at the dissipative scale, for which the Gaussian wavelet is no longer adapted.The results of the computations over 10 snapshots of velocity are shown in gure 4.10. The use of the mollied velocity gradient leads to smoother results

Figure 4 .

 4 Figure 4.10: (a) Joint PDF of Q and R. (b) Conditional average of the local energy transfer D I at given Q and R. (c) Conditional average of the local Hölder exponents at given Q and R. For each gures, three Vieillefosse lines are drawn in dashed lines.

Figure 4 . 11 :

 411 Figure 4.11: Time evolution of the energy dissipation rate as a function of time for the Run II. The rst vertical bar materializes the time at which the rst snapshot was extracted. The second vertical bar materializes the start of the time resolved dataset.

Figure 4 .

 4 Figure 4.12: Main gure: Multifractal spectra for the velocity elds from simulation in the dissipative regime. The red (resp. green, blue) curve corresponds to the global (resp. longitudinal, perpendicular) increments. The black curve is the result using the velocity increments based on singular value decomposition. The error bars correspond to a shift of the tting range by 10% for the power laws.

gure 4 . 14 .Figure 4 .

 4144 Figure 4.13: (a) Multifractal spectrum tted around ≈ 1.4λ. The spectrum corresponding to the global increments is shown in red, while the spectrum obtained using the unmodied WTMM method is shown in black for comparison. The vertical line materializes h = 1/3. The error bars correspond to a shift of the tting range by 5% for the power laws. (b) Box counting dimensions evaluated at the same scale. (c) Matching function resulting from the combination of (a) and (b).

Figure 4 .

 4 Figure 4.14: (a) Values of the threshold T p ( ) for the global increments as a function of scale for several p computed over Run II. The thresholds are renormalized by T 0 ( ). (b) Reference thresholds T 0 ( ). Fits corresponding to the inertial range, associated to ≈ 1.4λ are given in dotted lines while t in the dissipative range associated to = 1.8η are given in dashed lines.

Figure 4 . 15 :

 415 Figure 4.15: Matching function linking the velocity increment to a value of local Hölder exponent at the scale = 1.8η. The red curve and dots corresponds to the result of the calibration. The dotted blue line is the extrapolation used in case the velocity increment is higher than the 0.2.

( 4 . 2 )

 42 . We compute as before the joint PDF of the invariants Q and R and the conditional average for both the local Hölder exponents h and the local energy transfer D I . The results are shown in gure 4.18. The conclusions previously made in the inertial case are still applicable: We recover the common droplet shape for the joint PDF of Q and R in gure 4.18a. With the exception of a single excursion in the "vortex stretching" region, the strong positive Duchon Robert

Figure 4 . 16 :

 416 Figure 4.16: Joint PDF of h and the DuchonRobert energy transfer term D I over 100 elds from Run II. In order to enlighten the correlation of those terms independently of the number of occurrence, the PDF is rescaled such that every vertical line is a conditional PDF of D I at given h.

Figure 4 . 17 :

 417 Figure 4.17: Detail of the joint PDF of h and the DuchonRobert energy transfer term D I over 100 elds from run II corresponding to very low Hölder exponents. The PDF is rescaled such that every vertical line is a conditional PDF of D I at given h.

Figure 4 .

 4 Figure 4.18: (a) Joint PDF of Q and R. (b) Conditional average of the local energy transfer D I at given Q and R. (c) Conditional average of the local Hölder exponents at given Q and R. For each gures, three Vieillefosse lines are drawn in dashed lines.

Figure 4 . 19 :

 419 Figure 4.19: Histogram of the number of distinct events detected in the 40 velocity elds uncorrelated in time.

Figure 4 . 20 :

 420 Figure 4.20: Scatterplot of the mollied velocity gradient tensor invariants Q and R at the location of minimal Hölder exponent for the 213 extreme events extracted at the dissipative scale. Three Vieillefossse lines at 27R 2 + 4Q 3 = -62500, 0 and 62500 are drawn in dashed lines.

Figure 4 . 21 :

 421 Figure 4.21: Visualization of the 2nd strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 30% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 2 , 3 . (m): Isosurfaces of local energy dissipation at D ν = 5 , 6.5 .

Figure 4 . 22 :

 422 Figure 4.22: Visualization of the 15th strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 30% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 4 , 5 . (m): Isosurface of local energy dissipation at D ν = 5 , 6.5 .

Figure 4 . 23 :

 423 Figure 4.23: Visualization of the 30th strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 40% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 2 , 3 . (m): Isosurfaces of local energy dissipation at D ν = 4 , 5 .

Figure 4 .

 4 Figure 4.24: (a): Streamlines of velocity. The axial velocity is oriented away from the x = 0 plane. The red (resp. blue) isosurface corresponds to an helicity of H = 2.5 (resp. H = -2.5). (b): The magenta isosurface corresponds to the local Hölder exponent h = 1.35. The red isosurface corresponds to D I = 2 . The blue isosurface corresponds to D ν = 4 .

Figure 4 .

 4 Figure 4.25: (a) Fit of the component of the vorticity ω x in the x = 0 plane for the averaged event. The vorticity along the z (resp. y) direction is represented with blue circles (resp. red crosses) while the Burgers t is shown as a black continuous curve. (b) Colormap of the component of the vorticity ω x in the x = 0 plane for the averaged event.

Figure 4 . 26 :Figure 4 . 27 :Figure 4 . 28 :

 426427428 Figure 4.26: Time evolution of the local Hölder exponents (in magenta, left axis), the vorticity ωt ν (in green, middle axis), the local energy transfer D I / (in red, rightmost axis) and the energy dissipation D ν / (in blue, rightmost axis). The vertical lines materializes the ve snapshots represented in gures 4.27 and 4.28.

Figure 5 . 1 :

 51 Figure 5.1: Schematic representation in 2D of the particles being advected by the large scale velocity.

Figure 6 . 1 :

 61 Figure 6.1: (a) Test of universality of equation (3.49). The data points from run I are marked by triangles and the data points from run II are marked by circles. The orders of p are color coded. (b) Comparison of the coecients c p for the two runs. Blue circles correspond to Run I and red triangle to Run II.

Figure A. 1 :

 1 Figure A.1: Visualization of the 1st strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 30% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 2 , 3 . (m): Isosurfaces of local energy dissipation at D ν = 4 , 5 .

Figure A. 2 :

 2 Figure A.2: Visualization of the 11th strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 30% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 2 , 3 . (m): Isosurfaces of local energy dissipation at D ν = 4 , 5 .

Figure A. 3 :

 3 Figure A.3: Visualization of the 21st strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 30% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 2 , 3 . (m): Isosurfaces of local energy dissipation at D ν = 4 , 5 .

Figure A. 4 :

 4 Figure A.4: Visualization of the 31st strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 30% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 2 , 3 . (m): Isosurfaces of local energy dissipation at D ν = 4 , 5 .

Figure A. 5 :

 5 Figure A.5: Visualization of the 41st strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 30% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 2 , 3 . (m): Isosurfaces of local energy dissipation at D ν = 4 , 5 .

Figure A. 6 :

 6 Figure A.6: Visualization of the 51st strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 30% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 2 , 3 . (m): Isosurfaces of local energy dissipation at D ν = 4 , 5 .

Figure A. 7 :

 7 Figure A.7: Visualization of the 61st strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 30% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 2 , 3 . (m): Isosurfaces of local energy dissipation at D ν = 4 , 5 .

Figure A. 8 :

 8 Figure A.8: Visualization of the 71st strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 30% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 2 , 3 . (m): Isosurfaces of local energy dissipation at D ν = 4 , 5 .

Figure A. 9 :

 9 Figure A.9: Visualization of the 81st strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 30% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 2 , 3 . (m): Isosurfaces of local energy dissipation at D ν = 4 , 5 .

Figure A. 10 :

 10 Figure A.10: Visualization of the 91st strongest event detected. (a-c): Slices of velocity (the quiver plot codes for the vector in plane and the color for the third component.) (d-f ): Slices of vorticity (same as above) (h): Scatter plot of D I at coordinates (R,Q). (i): Scatter plot of h at coordinates (R,Q). (j): Streamlines of velocity. (k): Isosurfaces of vorticity at 30% and 60% of ωmax. (l): Isosurfaces of local Hölder at h = 1.3, 1.35. (m): Isosurfaces of local energy transfer at D I = 2 , 3 . (m): Isosurfaces of local energy dissipation at D ν = 4 , 5 .

FIG. 1 .

 1 FIG. 1. Main figure:Exponent τ as a function of the parameter q for the global (red circles), longitudinal (green squares), perpendicular (blue triangles) components as well as for the native method (black diamonds) applied to 100 fields of fBm in 2 dimensions with 3 components of velocity at H = 1/3. The 3 components are generated independently before enforcing the divergence free condition. The theoretical curve from equation (26) is materialized by the solid magenta line. Inset figure: Corresponding parametric plot of D(q) with respect to h(q). The dashed line materializes the theoretical value of H = 1/3.

FIG. 2 .

 2 FIG.2. Box counting dimension of the isocontours of the global (solid red), longitudinal (dashed green) and perpendicular (dotted blue) components of the wavelet based velocity increments for = 2.8 × 10 -3 for the fBm.

8 FIG. 3 .FIG. 4 .

 834 FIG.3. Values of the threshold Sp( ) for the global increments as a function of scale for several p for the fBm. The values of p increase from bottom to top.

FIG. 5 .

 5 FIG. 5.Pseudo "multifractal spectrum" reconstructed from the boxcounting dimensions of isovalues of velocity increments at scale = 8.2 × 10 -4 and the scaling laws over the Sp( ). The red circles (resp. green squares, blue triangles) correspond to the spectrum reconstructed using the global (resp. longitudinal, perpendicular) increments. The solid vertical line materializes the theoretical h = 1/3.

3 FIG. 6 . 3 FIG. 7 .

 3637 FIG.6. Energy spectra for the two simulations. The solid blue (resp. dashed red) curve stands for the simulation in the inertial (resp. dissipative) range. The black line materializes the k -5/3 slope for the inertial range. The horizontal lines of the same color as the curves materialize the corresponding fitting range used when computing power laws. The vertical lines materialize the k λ corresponding to each simulation.

3 FIG. 8 .

 38 FIG. 8. Main figure:Multifractal spectra for the velocity fields from simulation in the dissipative regime. The solid red (resp. dashed green, dotted blue) curve corresponds to the global (resp. longitudinal, perpendicular) increments. The dash-dotted black curve corresponds to the native method. The error bars correspond to a shift of the fitting range by 10% for the power laws. Inset: Exponents τ as a function of q for the global (red circles), longitudinal (green squares)and perpendicular (blue triangle) increments, as well as the native method (black diamonds).

FIG. 10 .FIG. 11 .

 1011 FIG.10. Values of the threshold Sp( ) for the global increments as a function of scale for several p for the simulation in the inertial range.

  FIG.12. Pseudo spectra reconstituted using the coefficient cp computed at scale ≈ 0.24 and the boxcounting dimensions of the velocity increments computed at scales ≈ 0.12 (dashed green curve), ≈ 0.24 (solid red curve) and ≈ 0.46 (dotted blue curve). The scale are also materialized by vertical lines of the same color on the Figure10.

FIG. 13

 13 FIG. 13. (a): Arrows stand for in plane velocity; color plot for the singularity exponent h(x). The gray areas corresponds to increments outside of the domain of the mapping function. (b): Arrows stand for in plane velocity; color plot for the Duchon-Robert term D I (x). (c): 3D representation of the same event. The green surface is an isosurface of h(x) at value h = 0.14, the red surface is an isosurface of D I (x) at value D I = 0.5 . All computations are done for a value of ≈ 0.24.

  We can assume that |Ω T | = 1. If V r,p = 1, the statement is trivial. Suppose V r,p < 1. Note that the function µ({|f | α}) is continuous from the left, and at a point of a jump the size of the jump is exactly µ({|f | = α}). Hence, there exists an α 0 and a set B ⊂ {|f | = α} such that A = {|f | > α} ∪ B has measure exactly V r,p . By Chebyshev, V r,p 1 α r A |f | r dµ.

A

  r,p = {|f | c r,p s r,p }. (B6) Directly by Chebyshev, we readily obtain the bound |A r,p | 1 c p r,p s p r,p |f | p |Ω T | = 1 c p r,p V r,p |Ω T |, and |f | p χ ΩT \Ar,p c p-r r,p |f | p |f | r |f | r χ ΩT \Ar,p c p-r r,p |f | p , which implies (1c p-r r,p ) |f | p |f | p χ Ar,p . (B7) Note again that as r → p, the information about concentration of f gets lost. However, the corresponding threshold amplitudes converge to a non-trivial value, namely, lim r→p s r,p = s p = exp |f | p ln |f | |f | p . (B8) Let us call it p-amplitude. One can easily restore s r,p from s τ by the formula s r,p = exp 1 pr p r ln s τ dτ .

A-5 3 /

 3 few years later, Frisch and Vergassola 4 claimed that the universality of the energy spectrum can be recovered, if the fluctuations of the dissipative length scale are taken into account by introducing a new non-dimensionalisation procedure. The new prediction is that log E(k)ǫ -2 3 η log(Re) should be a universal function of log(kη)/ log(Re), where Re is the Reynolds number. This claim was examined by Gagne et al., later using data from the Modane wind tunnel experiments 1 . They

FIG. 1 :× 10 4 ⋆

 14 FIG. 1: Von Krmn swirling flow generator. (a) normal view, bottom (b) and top (c) impellers rotating -both seen from the center of the cylinder, and (d) sketch with the relevant measures. A device not shown here maintains the temperature constant during the experiment. Both impellers are counter-rotating.

, ( 5 )

 5 This scaling is typical of a global scale symmetry solutions, and was criticized by Landau, who considered it incompatible with observed large fluctuations of the local energy dissipation. Kolmogorov then built a second theory (K62), in which fluctuations of energy dissipation were assumed to follow a log-normal statistics, and taken into account via an intermittency exponent µ and a new length scale L, thereby breaking the global scale invariance. The resulting velocity structure functions then follow the new scaling:S p (ℓ) ∼ (ǫℓ) p/3 ℓ L µp(3-p)which implies a new kind of universality involving the relative structure functions Sp as:

  τ (p) is the Legendre transform of the rate function C(h + 1/3), i.e. τ (p) = min h (p(h -1/3) + C(h)), and equivalently, that C(h) is the Legendre transform of τ (p). Because of this, it is necessarily convex. The set of points where it satisfies C(h) ≤ d, represents the set of admissible or observable h, is therefore necessarily an interval, bounded by -1 ≤ h min andh max ≤ 1.As noted by4 , the scaling exponent ζ(p) = p/3 + τ (p) defined via Equation (11) is only constant in a range of scale where ℓ > η h for any h ∈ [h min , h max ]. For small enough ℓ, this condition is not met anymore, since as soon as ℓ < η h , all velocity fields corresponding to h are "regularized", and do not contribute anymore to intermittency since they scale like ℓ. This results in a slow dependence of ζ(p) with respect to the scale, which is obtained via the corrected formula:S p = (ǫℓ) p/3 η h ≤ℓ ℓ L p(h-1/3)+C(h) dh ∼ (ǫℓ) p/3 ℓ L τ (p,ℓ).

  is simply: h(ℓ) = -1 + log(Re)/ log(L/η). We note θ = log(L/ℓ)/ log(Re). We can now rewrite equation -1/3)+C(h) dh ∼ exp (-θτ (p, θ) log(Re)) ,

γ(Re) log S p A p u p K=

 p G (p, γ(Re) log(ℓK 0 /η)) ,(15) where A p and K 0 are universal constant and β and G are general functions, F being linear in the inertial range, G(p, x) ∼ τ (p)x. The validity of this universal scaling was checked by Gagne and Castaing 1 on data obtained from the velocity fields measured in a jet using hot wire anenometry.They found good collapse of the structure functions at different Taylor Reynolds R λ , provided γ(Re) is constant at low Reynolds numbers and follows a law of the type:γ(Re) ∼ γ 0 / log(R λ /R * ),where R * is a constant, whenever R λ > 400. Since we have R λ ∼ Re 1/2 and (L 0 /η) ∼ Re 3/4 , we can rewrite equation (15) as:β(Re)log( Sp /S 0p ) log(L 0 /η) = G p, β(Re) log(ℓ/η) log(L 0 /η) ,

FIG. 2 :

 2 FIG. 2: Test of K41 universality Eq. (4). a) Numerical data b) Experimental data. The structure functions have been shifted by arbitrary factors for clarity and are coded by color: p = 1: blue symbols; p = 2: orange symbols; p = 3: yellow symbols; p = 4: magenta symbols; p = 5: green symbols; p = 6: light blue symbols; p = 7: red symbols; p = 8: blue symbols; p = 9: orange symbols. For K41 universality to hold, all the function should be constant, for a given p.

FIG. 3 :

 3 FIG. 3: Test of K62 universality Eq. (7). a) Numerical data b) Experimental data. The structure functions have been shifted by arbitrary factors for clarity and are coded by color: p = 1: blue symbols; p = 2: orange symbols; p = 3: yellow symbols; p = 4: magenta symbols; p = 5: green symbols; p = 6: light blue symbols; p = 7: red symbols; p = 8: blue symbols; p = 9: orange symbols. The dashed lines are power laws with exponents τ (p) = ζ(p)ζ(3)p/3, with ζ(p) shown in figure 4-a.

FIG. 4 :

 4 FIG. 4: . a) Determination of ζ(3) by best collapse using both DNS (open symbols) and experiments (filled symbols). The black dashed line is ℓ 0.8 . b) Scaling exponents ζ(p) of the wavelet structure functions of δW as a function of the order, from Table III, for DNS (blue circle) and experiments (red square) . The red dotted line is the function min h (hp + C(h)) with C(h) given by C(h) = (ha) 2 /2b, with a = 0.35 and b = 0.045. The black stars correspond to ζ SAW (p)/ζ SAW (3) (see Table III), while the black triangle correspond to ζ EXP (p)/ζ EXP (3).

FIG. 5 :

 5 FIG. 5: Test of general universality equation (20) using both DNS (open symbols) and experiments (filled symbols). The functions are coded by color. a) p = 1: blue symbols; p = 2: orange symbols; p = 4: magenta symbols; p = 5: green symbols; b) p = 6: light blue symbols; p = 7: red symbols; p = 8: blue symbols; p = 9: orange symbols. The functions have been shifted by arbitrary factors for clarity. The dashed lines are power laws with exponents τ (p) = ζ(p)ζ(3)p/3, with ζ(p) shown in figure 4-a.

3 FIG. 6 9 ζ

 369 FIG. 6: a) Variation of 1/β(Re) versus log(R λ ) in experiments (red square) and DNS (blue circle) when using the DNS at R λ = 138 as the reference case. We plotted in black the values found by Gagne and Castaing in 1 shifted by an arbitrary factor to coincide the values at large Reynolds. The black dashed line is (4/3) log(R λ /5). b) Multifractal spectrum C(h) for the experiments. The spectrum has been obtained by taking inverse Legendre transform of the scaling exponents ζ(p) shown in figure 4-b. The dotted line is a parabolic fit C(h) = (ha) 2 /2b with a = 0.35 and b = 0.045.

  . The exponents τ EXP (p)(red square) and τ DNS (blue circle) have been computed through a least square algorithm upon τ (p), minimizing the scatter of the rescaled structure functions log L 0 η τ (p) Sp ) with respect to the line (ℓ/η) τ (p) . The corresponding ζ(p) were inferred using the formula ζ(p) = τ (p) + ζ(3)p/3, where ζ(3) is computed in figure 4-a. 6-b for both the DNS and the experiment.

3 =

 3 where S(E) is the large deviation function of log(µ ℓ ) and has the meaning of an energy while log(ℓ/η) has the meaning of a volume, and log(µ ℓ )/ log(ℓ/η) is an energy density. Because of the definition of µ ℓ , it is easy to see that S is connected to C, the large deviation function of |δW ℓ . In fact, since in the inertial range where|δW ℓ | 3 ∼ ℓ ζ(3) , we have S(E) = C(3h-ζ(3)). By definition, e p log(µ ℓ ) ,

FIG. 8 :

 8 FIG.8: Same as figure7for R λ = 90 (line), R λ = 138 (dotted line). Note the inflexion point appearing in the curves.

Table 4 .

 4 1: Parameters used in the simulations with grid resolution 768 3 . ν is the kinematic viscosity, L is the integral scale, λ is the Taylor scale, η is the Kolmogorov scale, u 2 1/2 is the rms velocity, is the energy dissipation rate, ηk max characterizes the resolution (ηk max > π is well resolved), Re L is the Reynolds number based on the integral scale, and Re λ is the Reynolds number based on Taylor scale λ.

	Run	ν	L	λ	η	u 2 1/2	ηk max Re L Re λ
	I II	7.5 × 10 -4 0.79 0.19 0.0083 5 × 10 -3 0.94 0.48 0.034	0.54 0.55	0.089 2.1 0.097 8.5	570 104	140 53

30 eddy turnover times. As a result, the elds are uncorrelated, which ensures better statistical convergence.

TABLE I .

 I yielded a most probable Hölder Parameters used in the simulations. N is the linear grid resolution, ν the kinematic viscosity, L is the integral scale, λ is the Taylor scale, η is the Kolmogorov scale, u 2 1/2 is the rms velocity, is the energy dissipation rate, ηkmax characterizes the resolution (ηkmax > π is well resolved), Re is the Reynolds number based on the integral scale, and Re λ is the Taylor based Reynolds number.

	Run	N	ν	L	λ	η	u 2 1/2		ηkmax	Re	Re λ
	I II	768 3 768 3	7.5 × 10 -4 5 × 10 -3	0.79 0.94	0.19 0.48	0.0083 0.034	0.54 0.55	0.089 0.097	2.1 8.5	570 104	140 53

  FIG.9. Box counting dimension of the isosurfaces of the different components of the wavelet based velocity increments for ≈ 0.24 for the simulation in the inertial range. The color code for the different components is the same as for the previous figures.

		3.0
		2.8
		2.6
	D bc (T )	2.2 2.4
		2.0
		1.8
		1.6
		0.0	0.1	0.2	0.3	0.4	0.5	0.6
					T

  with 0 < c p < 1 and s p as above captures an appreciable concentration of the renormalized probability density F p = |f | p |f | p , provided U 0 as chosen appropriately. So, let us define as before A p = {|f | c p s p }.

	Let us consider the entropy concentration on the comple-
	ment of A p :
	|f | p ln	|f | U 0	χ ΩT \Ap	ln c p |f | p χ ΩT \Ap
						+	|f | p ln |f | U0 |f | p	|f | p χ ΩT \Ap (B12)
	Let us choose U 0 = |f | p 1/p , and denote I = |f | p ln |f | U0 . Then
					p |f | p I = F p ln F p ,
	which by the classical Csiszar-Kullback bounded from
	below by F p -1 2 1 , i.e. non-negative. This allows to continue the line of (B12):
	|f | p ln	|f | U 0	χ ΩT \Ap	ln c p |f | p χ ΩT \Ap + |f | p ln	|f | U 0	.
	(B13) |f | p and rewriting Multiplying the above by the factor of p in terms of F p we obtain
		F p ln F p χ Ap	F p χ ΩT \Ap ln	1 p c p	.	(B14)
	Note that ln 1 c p
	(B11)

  Table I summarizes the different parameters; R λ = λu rms /ν is the Reynolds number based on the Taylor length scale λ =

	u 2

TABLE II :

 II 6 , whereas the run with R λ = R λ η k max η N x × N y ×N z ℓ min /η k f Samples Symbol Parameters for the DNS. R λ is the Reynolds based on the Taylor micro-scale. η is the Kolmogorov length. The third column gives resolution of the simulation through k max η,where k max = N/3 is the maximum wavenumber. The fourth column gives the grid size; notice that the length of the box is 2π. Here, ℓ min is the smallest scale available for the calculations of the wavelets. kf is the forcing scale. The Sample columns gives the number of points (frames × gridsize) over which are calculated the statistics.

	25 0.079 3.35	128 3	0.635 1 5000	⋆
	56 0.034 6.42	256 3	0.31 1 105 000 △
	80 0.020 1.68	256 3	1.22 1 270 000
	90 0.017 5.70	1024 3	0.36 1 10 000 ♦
	138 0.009 1.55	512 3	1.37 1 12 000 •

TABLE III :

 III Scaling exponents τ (p) and ζ(p) found by the collapse method based on K62 universality for experimental data (subscript EXP) or numerical data (subscript DNS). The subscript SAW refers to the values obtained by 5 on the same set of experimental data, using velocity increments and Extended Self-Similarity technique11

TABLE IV :

 IV Summary of the analogy between the multifractal formalism of turbulence and thermodynamics.The thermodynamic analogy is thus meaningful and is summarized in TableIV. It can be used to derive interesting prospects.

p > 0. Thus, (B14) expresses quantitatively a concentration proportion of the entropy relative to the complementary probability of the event A p .

Reproduction of an article in Entropy (2019).

P. Debue, V. Shukla, D. Kuzzay, D. Faranda, E.-W. Saw, F. Daviaud, and B. Dubrulle,

Remerciements Acknowledgments. FN has been funded by École Normale Supérieure and was granted access to the HPC resources of IDRIS under the allocation 2A10096 made by GENCI (Grand Équipement National de Calcul Intensif). This work was funded by the ANR EXPLOIT, grant agreement no. ANR-16-CE06-0006-01. Research was supported in part by NSF grant DMS 1515705, Simons Foundation, and the College of LAS at UIC.

Ce travail a été nancé par le biais d'une Allocation Spécique Normalienne,

Appendix A

Visualization of events

This appendix is dedicated to the visualization of several more singular events extracted from Run II. The selection is made to cover a large range of events, not restriced to the few most singular ones.

Local estimates of Hölder exponents in turbulent vector elds

Reproduction of an article in Physical Review E (2019).

Appendix A: General considerations about velocity increments

Consider a 3D velocity field u. Without loss of generality, it can be decomposed using Helmoltz decomposition as:

By construction, ∇ × u = -∆Q and ∇ • u = ∆φ. So, if the fluid is incompressible (which we do not assume for the present time), ∆φ = 0. We can further define velocity increments of such field over a scale as:

that can be further decomposed into longitudinal δ L u( ) and perpendicular δ P u( ) velocity increments as:

With this definition, δ L u(r) is a scalar and δ P u(r) is a vector orthogonal to both δu and r. The definition of the perpendicular velocity increment is similar to the more commonly used transverse increment, as they have the same norm but different orientation.

We consider now the angle average of r δ L u(r) and r δ P u(r) over a sphere of radius , defined through convolution with a Gaussian function Ψ(x) = exp(-x 2 /2)/N , where N is chosen so that the Gaussian is of norm unity:

Using properties of the Gaussian, it can be easily checked that:

Using this property and the Helmoltz decomposition, we then see that after an integration by part that:

Rephrasing this, we see that the angle averaged longitudinal velocity increment is connected with the potential part of the Helmoltz decomposition (or the divergence of the velocity field), while the angle averaged transverse velocity increment is connected with the rotational part of the Helmoltz decomposition (or the vorticity). The connection is through a wavelet transform, using a Gaussian wavelet. That way, studying the scaling properties of either δ L u(r) r ang, or δ P u(r) r ang, provides information about the scaling of the potential or rotational part of the velocity. Such decomposition has already been used by [14] and [15] to study structure functions in experimental turbulence.

Appendix B: Volumetric concentration and amplitude factors

In this section we introduce several scalar factors associate to a general vector or scalar field f , which capture a measure of concentration and size of f . We use the following definition of the space-time average

where Ω T = Ω × [0, T ], and Ω denotes a fluid domain.

We also denote dµ = 1 |ΩT | dx dt, the averaging measure.

Concentration factors

We now present formalities of the active volume/region theory. Let f ∈ L ∞ (Ω T ) be an arbitrary field. Let 1 p, r ∞, r = p. We define the (r, p)-concentration factor of f as follows

Let us list a few elementary properties: V r,p is adimensional, V r,p = V p,r , and by Hölder, V r,p 1. For any triple p 1 < p 2 < p 3 we have by interpolation,

The idea is that V r,p measures a proportion of the volume of Ω T where much of the L p -weight of f concentrates provided p > r. More precisely, we have the following lemma. We investigate universality of the Eulerian velocity structure functions using velocity fields obtained from the stereoscopic particle image velocimetry (SPIV) technique in experiments and the direct numerical simulations (DNS) of the Navier-Stokes equations. We show that the numerical and experimental velocity structure functions up to order 9 follow a loguniversality 1 ; we find that they collapse on a universal curve, if we use units that include logarithmic dependence on the Reynolds number. We then investigate the meaning and consequences of such log-universality, and show that it is connected with the properties of a "multifractal free energy", based on an analogy between multifractal and themodynamics.

We show that in such a framework, the existence of a fluctuating dissipation scale is associated with a phase transition describing the relaminarisation of rough velocity fields with different Hölder exponents. Such a phase transition has been already observed using the Lagrangian velocity structure functions, but was so far believed to be out of reach for the Eulerian data. ii) Using the general universality, we may also get τ p,univ by a linear regression on the collapse curve.

Note that since the data are collapsed, this provides a very good estimates of this quantity, with the lowest possible noise. In practice, we observe no significant differences with the two estimates; therefore, we only report the values obtained by following the first method.

The corresponding values are plotted in figure 4

Practically, this amount to use the following formula:

To estimate C, we thus first perform a polynomial interpolation of order 4 on ζ(p), then derivate the polynom to estimate d ζ(p) d p , thus get C through equation ( 22). The result is provided in figure

B. Multifractal pressure and phase transition

Given our free energy, F = log( S3p ), we can also compute the quantity conjugate to the volume, i.e. the multifractal pressure as: P = ∂F/∂V . In the inertial range, where Sp ∼ ℓ τ (p) , we thus get P = τ (p), which only depends on the temperature. Outside the inertial range, P has the meaning of a local scaling exponents also depends upon the scale, i.e., on the volume V and on N (Reynolds number). Using our universal functions derived in figure 5, we can then compute empirically the multifractal pressure P and see how it varies as a function of T , V and N . It is provided in figure 7 for R λ = 25 and R λ = 56, and in figure 8 for R λ = 90 and R λ = 138. We see that at low Reynolds number, the pressure decreases monotonically from the dissipative range, reaches a lowest points and then increases towards the largest scale. There is no clear flat plateau that would correspond to an "inertial" range. In contrast, at higher Reynolds number, a plateau appears for p = 1 to p = 4 when going towards the largest scale, the value of the plateau corresponding to τ DNS . The plateau transforms into an inflection point for p ≥ 5 making the derivative ∂P/∂V change sign. This is reminiscent of a phase transition occurring in the inertial range, with coexistence of two phases: one "laminar" and one "turbulent". We interpret such a phase transition as the result of the coexistence of region of flows with different Hölder exponents, with areas where the flow has been relaminarized due to the action of viscosity, because of the random character of the dissipative scale (see below).

VI. CONCLUSION

We have shown that a deep analogy exists between multifractal and classical thermodynamics. In this framework, one can derive from the usual velocity structure function an effective free energy that respects the classical extensivity properties, provided one uses a number of degrees of freedom (given by N = 1/β(Re)) that scales like log(R λ ). This number is much smaller than the classical N ∼ Re 9/4 that is associated with the number of nodes needed to discretize the Navier-Stokes equation down to the Kolmogorov scale. It would be interesting to see whether this number is also associated with the dimension of a suitable "attractor of turbulence". Using the analogy, we also found the "multifractal" equation of state of turbulence, by computing the multifractal pressure P = ∂F/∂V . We found that for large enough R λ and p (the temperature), the system obeys a phase transition, with coexistence of phase like in the vapor-liquid transition. We interpret this phase transition as the result of the coexistence of region of flows with different Hölder exponents,