Jean-Paul Caltagirone
email: calta@ipb.fr

Turbulence initiation in a Taylor-Green flow by vortex bursting in an edge-based formulation

Keywords: Turbulence cascade, Vortex bursting, Discrete Mechanics, Conservation of acceleration, Helmholtz-Hodge Decomposition, Dynamic entanglement

The development of turbulent structures and their decay are analyzed from a physical model based on discrete mechanics. The associated acceleration-preserving equation of motion establishes that the acceleration of the material medium is equal to the two accelerations applied to it, the first expressing the compressive effects and the second the viscous effects. The equation is thus presented as a Helmholtz-Hodge decomposition in two terms, one with divergence-free and the other with curl-free. The case treated in direct simulation, the Taylor-Green vortex, is widely used to test the accuracy and the performance of high-order methods but especially to analyze turbulent flows with high Reynolds numbers solved by the Navier-Stokes equation. While the results of the Navier-Stokes equations and discrete mechanics are strictly the same for many classical solutions, they differ significantly for the Taylor-Green vortex flow. The energy decreases monotonically with the Navier-Stokes equation while it increases in the inertial region until the onset of vortex bursting which triggers the energy cascade and viscous dissipation process. A careful verification of the purely numerical aspects allows to rule out possible artifacts and to conclude that the differences are due to the physical models implemented. An in-depth physical analysis allows to understand the origin of the differences between the two physical models.

Introduction

The direct simulation of turbulence is closely linked to the Navier-Stokes equation, the main way of the existing continuous medium framework to approach the behavior of vortex structures from their creation to their dissipation [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF][START_REF] Taylor | Mechanism of the production of small eddies from large ones[END_REF][START_REF] Frish | Turbulence, the legacy of A.N. Kolmogorov[END_REF]. Since the middle of the 18th century, the equation derived from Euler's equation associates the effects of inertia and compression with the notion of viscosity and dissipation.

Emblematic cases of turbulence such as Homogeneous Isotropic Turbulence [START_REF] Sagaut | homogeneous turbulence dynamics[END_REF], the turbulence in fully developed channel flow [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF], or the decay of the turbulence in Taylor-Green vortex allow to study the physical characters of this one. This last case was approached several decades ago using pseudo-spectral methods of very high precision [START_REF] Orszag | Analytical theories of turbulence[END_REF][START_REF] Gottlieb | Numerical Analysis of Spectral Methods: Theory and Applications[END_REF], it has since become a relevant test case to verify the accuracy and robustness of numerical methods. This will be the case chosen to evaluate the physical model associated with the discrete mechanics [START_REF] Caltagirone | Discrete Mechanics, concepts and applications[END_REF].

The discrete formulation is based on a geometrical approach of space where the notion of continuous medium is abandoned; the Galilean frame of reference of classical mechanics is replaced by a local frame of reference where the long distance interactions are of cause and effect. In the same way, time is replaced by a discrete notion where the knowledge of the only solution at a given time t o allows to predict the one at a time t o + dt where dt is the time lapse between two mechanical equilibria. This physical model corresponds to a continuous memory view where the history of motion is stored in potentials. The notion of mass is put aside in favor of a description based on the total energy represented by the discrete equation of motion [START_REF] Caltagirone | An alternative to the concept of continuous medium[END_REF]. An essential difference is the dynamic entanglement of compression and rotation effects suggested by J.C. Maxwell [10] to federate magnetism and electrodynamics within what has been called electromagnetism.

However, this quite different approach has allowed, for fluid flows, to find exactly a large number of solutions of the Navier-Stokes equation. The analytical solutions or the result of simulations performed with both formulations are the same. Turbulence is of a different nature, it reveals different physical mechanisms, including the transfer of energy between large structures of the flow to smaller and sometimes the reverse. It is therefore an opportunity to test discrete mechanics to assess its relevance to a problem that remains one of the key elements of fluid physics.

After a brief presentation of the basics of discrete mechanics, direct simulations of the Taylor-Green vortex are presented in the form of time spectra of the kinetic energy, its derivative, the rotation energy, the compression energy, etc. The second part is devoted to a study of convergence in space of the turbulent statistics in order to remove any doubt on the quality of the numerical formulation associated with the discrete mechanics. Finally the differences between the results of the Navier-Stokes model and the discrete model are analyzed in depth. The dissimilarity of the equations on physical effects, inertia, rotational effects, dynamical entanglement allows to make sense of the differences observed on the results.

Discrete mechanics framework

Physical principles

The principles of physics and mechanics are not called into question; they will simply be revisited from the discrete angle chosen to derive an equation of motion which is presented as an alternative to the equations of the mechanics of continuous media.

One of the oldest principles, dating from the beginning of the 16th century and initiated in particular by Galileo, expresses the equivalence between inertial mass and gravitational mass; the acceleration of a body is then independent of its nature. This principle, since titled the Weak Equivalence Principle, was translated by Newton in the form of his second law, the fundamental law of dynamics, m γ = F [START_REF] Newton | Principes Mathématiques de la Philosophie Naturelle traduit en francais moderne d'après l'oeuvre de la marquise Du Châtelet sur les Principia[END_REF]. This logical interpretation in the historical context, where gravitation was the dominant theory of the preoccupations of the time, still exists today in many theories such as that of relativity, fluid flows, etc. Mass is attached to other quantities like energy (e = m c 2 ) or momentum q = m v. The Navier-Stokes equation is presented as a conservation of momentum. The interpretation given in discrete mechanics is different: if the mass does not intervene on the acceleration of a body subjected to an external acceleration, then why keep it? The postulate of discrete mechanics is to extend the equality of intrinsic and gravitational accelerations to any other acceleration. The law of discrete dynamics expresses the fact that the acceleration of a material medium or of a particle with or without mass is equal to the sum of the external accelerations; it is written:

γ = h ( 1
)
where h is the sum of all accelerations applied to the particle or to the material medium. This law is intrinsically valid independently of the velocity of the medium, whether the latter is in uniform rectilinear translational movement or at rest. The law [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF] expresses the conservation of acceleration on a rectilinear support. Acceleration is a quantity that can be measured anywhere, it is an absolute quantity, while velocity, force, energy, etc. are relative quantities, the knowledge of which in a previous state is required.

The Galileo principle of relativity, where all the laws of mechanics are identical in all the Galilean frames of reference, is interpreted in a different way since the notion of global frame of reference has been abandoned. The current velocity is here set from the expression v = v o + γ dt, where v o is the velocity at the initial time and dt the considered time increment. If the acceleration is zero, the velocity is kept on a rectilinear path. The absolute velocity of a body at uniform velocity is not required, it must naturally disappear from the equation of motion. The uniform rotational motion, which is a little different, also satisfies this principle of relativity under certain conditions [START_REF] Caltagirone | On Helmholtz-Hodge decomposition of inertia on a discrete local frame of reference[END_REF]. The discrete law of motion must then integrate the temporal process necessary to deduce the state of a system at time t from that at time t o ; the proposed model is continuous memory.

The notion of tensor initiated by mechanics for mechanics is also abandoned. It was linked to the nature of materials whose properties, depending on the directions of space, can be expressed simply by tensors. In discrete mechanics the notion of global frame of reference is replaced by a local frame of reference on a discrete geometric topology.

The principle of causality must of course be satisfied in the development of any physical model. In discrete mechanics this principle is satisfied if we consider that a piece of information can only be perceived over a distance dh, the discrete horizon, if dt is greater than dt = dh/c. In a system with several dimensions it is important to know how information is transmitted in a material medium, as this conditions the entanglement of the different mechanical effects.

Finally, the Helmholtz-Hodge decomposition, which is not really a principle, consists in separating the curl-free component and the divergence-free component of any vector. In mathematics, the Helmholtz-Hodge decomposition is useful for projecting a vector onto a space with zero divergence [START_REF] Bhatia | The Helmholtz-Hodge Decomposition -A Survey[END_REF][START_REF] Ranocha | Discrete Vector Calculus and Helmholtz Hodge Decomposition for Classical Finite Difference Summation by Parts Operators[END_REF]. In fluid mechanics, it easily ensures the incompressibility constraint [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF][START_REF] Angot | A fast vector penalty-projection method for incompressible non-homogeneous or multiphase Navier-Stokes problems[END_REF]. Numerous publications on the Helmholtz-Hodge decomposition [START_REF] Ahusborde | A primal formulation for the Helmholtz decomposition[END_REF][START_REF] Bhatia | The Helmholtz-Hodge Decomposition -A Survey[END_REF][START_REF] Angot | Fast discrete Helmholtz-Hodge decomposition in bounded domains[END_REF] show its theoretical interest and the potentialities of its applications. This decomposition of acceleration a priori is one of the cornerstones of discrete mechanics, it makes it possible to separate the terms of compression from the terms of rotation describing the motion of a material medium.

The physical parameters are essentially the longitudinal and transverse celerities, which must be located precisely within the discrete frame of reference. Indeed, the framework of classical mechanics is not adapted to a discrete vision based on these physical principles; in particular because of the concept of continuous medium where all the quantities are defined in a point. For this reason it is essential to define a discrete geometric topology from which the basic phenomena are described.

Primal and dual geometric structures

The derivation of the discrete equation of motion requires a complete geometric framework; the notions of continuous medium, derivation at a point, analysis and continuum mechanics are abandoned and replaced by two topologies called primal and dual, which can be inverted. The privileged directions of classical mechanics, like the notion of global frame of reference, are also abandoned. Figure [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF] shows schematically the two elementary geometric topologies; their orientation in space is not specified -one will simply define a local frame of reference by the orientations of the unit vectors t of the segment Γ and n orthogonal to the primal facet S.

The oriented contour Γ * of the primal facet is formed by the set of sides of the corresponding polygon; the primal S and dual ∆ facets are considered to be planar, even if this is not always the case in practice for general unstructured meshes. The contour of the dual facet is made up of all the oriented segments δ of the associated polygon. The vectors t and n are orthogonal by construction, t • n = 0. The segment Γ of length d = [a, b], whose ends are a and b, is chosen to represent the component of velocity v and the component of acceleration γ. The velocity vector and the acceleration vector in space will remain undefined; even if they can be reconstructed using their components, this will not be useful. Two functions will be used subsequently, φ a scalar defined at points a and b named scalar potential and ψ associated with the unit vector n named vector potential.

It is necessary to come back to the definition of the quantity v as well considered as a vector, a component of a vector, as a scalar on an oriented segment. This quantity is constant on Γ, it represents the average of the local velocity on the segment which varies from v a to v b . Assuming a linear variation, v is the arithmetic mean and, in general, it represents the mean value in the integral sense. The scalar |v| 2 = v • v is defined on the segment but also on each vertex |v a | 2 and at the center of the primal facets S. In addition, four discrete differential operators are introduced. The first is the operator ∇φ applied to the scalar potential; it differs from the classical gradient vector because it only applies to the component on the segment Γ, i.e. ∇φ = (φ b -φ a )/d. Thus ∇φ is both a scalar on the oriented segment and a vector. The second operator, called primal curl, is computed from the circulation of the vector v on the contour Γ * using Stokes' theorem, ∇×v; it will then be assigned to a dual segment δ and named ψ, one of the curl components orthogonal to n. Similarly, the dual curl operator ∇ d × ψ corresponds to a line integral on the dual contour, then assigned to the primal segment Γ. Finally, the divergence operator corresponds to the integral on the dual volume of all the fluxes having the same vertex, ∇ • v; the result is then assigned to the vertex a. In a continuous context, v is a polar vector and ψ an axial or pseudo-vector. The two primal and dual curl operators are thus differentiated by the function to which they apply. This discrete framework is completed by two notions, (i) the time lapse dt = t -t o between the current time and the time t o where the initial equilibrium is defined, and (ii) the discrete horizon dh = c dt, the length traveled by a wavefront of sound celerity c on a rectilinear trajectory.

The physical model has many properties, in particular the local and global orthogonality of the operators, ∇φ • ∇ d × ψ = 0, also verified by the discrete operators; the discrete model also mimics the properties of the continuum ∇ • ∇ d × ψ = 0 and ∇ × ∇φ = 0, whatever the polygonal or polyhedral geometric topologies, structured or unstructured [START_REF] Lipnikov | Mimetic finite difference method[END_REF][START_REF] Caltagirone | On primitive formulation in fluid mechanics and fluid-structure interaction with constant piecewise properties in velocity-potentials of acceleration[END_REF].

The notion of scalar potential is very general, φ can represent pressure, an electric potential, a gravitational potential, inertial, capillary, etc. It keeps its meaning in all media including vacuum; this is not the case of pressure which tends towards zero at the same time as the density. In the following section on fluid flows, the potential is equal to φ = p/ρ γ where γ is the ratio of the heats of mass.

The abandonment of the global reference frame of classical mechanics in favor of a local reference frame changes our perception of reality and of the physical model to represent it. The loss of the space vector abstraction is compensated by the definition of precise geometrical entities. In a way, a path only exists if it is drawn, the notion of free space vector is abandoned. Thus, all quantities, and in particular acceleration, are linked to a single segment Γ whose orientation remains undetermined. All interactions are defined by cause and effect.

Conservation of acceleration

The conservation of acceleration on a segment Γ, considered as a postulate of discrete mechanics, can be strongly related to the conservation of energy. Indeed, acceleration γ on segment Γ in figure [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF] is constant and its integration between a and b is none other than the variation of the total energy per unit of mass Φ between these two points:

Φ b -Φ a = b a γ • t dl (2)
Even if the velocity is not zero, the fact that it is also constant on the segment leads to zero energy variation according to the principle of relativity; all the terms of the particle derivative, the time derivative and the inertia are indeed equal to zero. The equivalence between mass and energy of the theory of relativity must be verified in this discrete context. It is therefore possible to impose the conservation of energy without having to additionally require the conservation of mass. The choice fixed in this context is to conserve the total energy and to eliminate the notion of mass, in accordance with the principle of discrete equivalence described above.

The acceleration vector, or its component γ, is decomposed into a curl-free part and another divergence-free part according to the Helmholtz-Hodge decomposition:

γ = -∇φ + ∇ d × ψ (3)
where φ is the scalar potential of the acceleration located at the points of the primal topology and ψ its vector potential assigned on the segment δ of normal unit n.

γ • t dl = - Γ ∇φ • t dl + Γ ∇ d × ψ • t dl (4)
where φ is the scalar potential of the acceleration located at the points of the primal topology and ψ its vector potential assigned on the segment δ of normal unit n.

The local (3) and integral (4) laws thus represent the conservation of energy and consequently the conservation of mass. The direct action represented by -∇φ cannot be modified by the induced action defined by ∇ d × ψ, or vice versa, because these two components are orthogonal. The entanglement resides in the fact that the intrinsic acceleration γ is necessary to allow the exchanges between the two orthogonal components. This entanglement does not depend on the spatial scale and can act at a great distance, but the effects depend on celerities c l and c t .

Several equations in physics present this dynamic entanglement, in solid or fluid mechanics or in electromagnetism. Some do not, for example Darcy's equation which derives only from a scalar potential φ D , such that v = -∇φ D . The generalized Darcy equation ∂v/∂t = -∇p -µ/K v becomes ∂v/∂t = -∇ (p + φ D ). The viscous effects act at a very short spatial scale of order of magnitude of √ K, the permeability of the porous medium.

Discrete equation of motion

The derivation of the discrete law of motion [START_REF] Caltagirone | Discrete Mechanics, concepts and applications[END_REF] serves to establish an equation to predict the solution of a problem at current time t from that at time t o , which represents the state of the physical system in mechanical equilibrium at an earlier instant. The discrete law defines a continuous memory model where the current mechanical equilibrium is calculated by a process of accumulation of compressive and shear stresses. The potentials at instant t o , φ o and ψ o , are energies per unit mass accumulated between the initial instant and instant t o ; they are called retarded potentials with reference to those of Liénard and Wichert [21]. The local equation of motion is then written:

dv dt = -∇ (φ o + dφ) + ∇ d × (ψ o + dψ) (5)
The increments of potentials dφ and dψ must be modeled with fundamental compression and shear experiments in mind. The variation of the scalar potential is related to the compression characterized by the divergence of the velocity and to the longitudinal celerity in the form dφ = -dt c 2 l ∇ • v. The variation of the vector potential depends on the curl of the velocity and the transverse celerity dψ = -dt c 2 t ∇ × v. The balance of the acceleration (4) on the segment Γ becomes:

Γ γ • t dl = - Γ ∇ φ o -dt c 2 l ∇ • v • t dl + Γ ∇ d × ψ o -dt c 2 t ∇ × v • t dl (6)
The retarded potentials reflect the storage of energy from the initial moment, or at least from a previous state of mechanical equilibrium, at moment t o . The history of compressive and shear energies is represented by these two retarded potentials:

φ o = - t o 0 c 2 l ∇ • v dτ ; ψ o = - t o 0 c 2 t ∇ × v dτ (7)
The equation of local motion can then be extracted by considering that distance d = [a, b] tends to zero. However, unlike the notion of continuous medium where all the quantities are reduced to a point, discrete mechanics interprets this step as a homothetic decrease towards zero of the distances, which preserves the angles between the segments of the geometric topologies. The law of motion and the temporal updates of potential reads:

           γ = -∇ φ o -c 2 l dt ∇ • v + ∇ d × ψ o -c 2 t dt ∇ × v + h s α l φ o -c 2 l dt ∇ • v -→ φ o α t ψ o -c 2 t dt ∇ × v -→ ψ o (8)
where h s represents the source terms, gravitational acceleration, capillary acceleration, etc. The intrinsic acceleration of the material medium or of the particle will be expressed from the partial derivative in time and the terms of inertia. The quantities α l and α t between 0 and 1 are attenuation factors for longitudinal and transverse waves; when α = 0, the waves persist indefinitely.

The symbol -→ corresponds to the temporal update of the quantity associated between t o and t; it is a discrete integration which traces the accumulation of mechanical stresses over time.

The material derivative dv/dt is modeled in a different way from those classically admitted in continuum mechanics; it was the subject of a specific presentation in [START_REF] Caltagirone | On Helmholtz-Hodge decomposition of inertia on a discrete local frame of reference[END_REF]:

γ = ∂v ∂t + ∇ Å 1 2 |v| 2 ã -∇ d × Å 1 2 |v| 2 n ã (9)
The inertial terms are described by a Helmholtz-Hodge decomposition of the inertial potential φ i = |v| 2 /2. The sum of these potentials with those of the acceleration will be called Bernoulli potentials, φ o B = φ o + φ i and ψ o B = ψ o + ψ i where ψ i = φ i n. Moreover the sum of the two terms of inertia ( 9) ∇φ i and ∇ d × ψ i can be interpreted as the curvature in space of the inertial potential [START_REF] Caltagirone | On Helmholtz-Hodge decomposition of inertia on a discrete local frame of reference[END_REF]. The intrinsic acceleration γ = ∂v/dt + ∇φ i -∇ d × ψ i introduces the notion of time thus closely coupling not only the inertial effects but also the compression and rotation effects in an entwined formulation.

Unlike longitudinal waves, transverse waves are polarizable in any planar surface orthogonal to the primal facet around the unit normal n. The shear in a given direction in the surface of the facet S induces a transverse wave of celerity c t in the orthogonal direction n. For example, for the plane startup Couette flow generated by a wall, the propagation of the transverse waves takes place in the direction orthogonal to the wall. In fluid, the depth of penetration of the waves is limited by viscous dissipation. It should be noted that for the Couette flow, the particular form of inertia in discrete mechanics leads, as in continuum mechanics, to zero inertia but only by compensation; the last two terms of the relation ( 9) are indeed different from zero. As the Reynolds number increases this flow becomes turbulent. In the general case of turbulent flows, the time constants of the vortex structures decrease and the interactions between the effects of pressure and shear are greater; these effects are closely intertwined with inertia.

The velocities c l and c t are intrinsic properties of the media considered (solid, fluid, vacuum) but they can of course vary according to other physical quantities such as temperature for example; they must simply be fixed in time and space. For fluids the product dt c 2 t is replaced by the kinematic viscosity ν for large time constants such that t > τ where τ ≈ 10 -11 s is the characteristic attenuation time of transverse waves. In the following the attenuation factor will be fixed at unity, α t = 1. The instantaneous vector potential ψ o = ν ∇ × v is not zero but the accumulation of rotational stresses is non-existent for Newtonian fluids.

The incompressibility constraint ∇ • v = 0 is never imposed in discrete mechanics; however, the approximation can be obtained for large values of the velocity. For values of order one of velocity and space we have ∇ • v ≈ c -2 l ; thus for values of c l larger than a few tens the flow can be considered as incompressible. The increase of c l towards 10 3 -10 8 ensures the incompressibility condition at machine precision. It should be noted that the Mach number M = v/c l is insufficient to decide whether a motion is compressible or not.

In fact it is the grouping dφ = dt c 2 l ∇ • v which defines the increase in compressive energy from one state of mechanical equilibrium to another and if c l increases the divergence decreases but dφ remains practically constant. Thus the variation of the compression energy with time is not negligible even for a flow considered incompressible. This is one of the cornerstones of discrete mechanics. The effects of compression in turbulent flows known as incompressible are too often neglected.

The particular symmetrical form of the equation of motion [START_REF] Caltagirone | Discrete Mechanics, concepts and applications[END_REF] in two Lagrangians induces according to Noether's theorem [START_REF] Kosmann-Schwarzbach | Noether Theorems. Invariance and Conservations Laws in the Twentieth Century[END_REF], the invariance of certain quantities, in particular the mechanical energy linked to temporal independence of physical law.

Discrete kinetic energy theorem

The kinetic energy theorem is not a new equation, it is the transformation of the discrete equation of motion [START_REF] Caltagirone | Discrete Mechanics, concepts and applications[END_REF] obtained by multiplying it by v. This is a scalar product of vectors having the same support Γ. Indeed ∇φ and ∇ d ×ψ are the two components of the acceleration γ carried by Γ. Considering the case of the fluid for which α l = 0, α t = 1 and ψ = -ν∇ × v, the kinetic energy theorem becomes:

1 2 d|v| 2 dt = -v • ∇φ o -v • ∇ d × (ν ∇ × v) (10)
Its integral form is obtained by integrating (10) on Γ:

1 2 Γ d|v| 2 dt dl = - Γ v • ∇φ o dl - Γ v • ∇ d × (ν ∇ × v) dl (11)
Since the quantity v is a component of the velocity vector remaining indeterminate and ∇φ o is the restriction on Γ of the classical gradient of a scalar, the vectors v and ∇φ are collinear and the scalar product can only be zero if one of the two vectors is equal to zero. The same is true for the second term on the right where ∇ d × (ν ∇ × v) is the result projected on Γ of the circulation of ψ along the dual contour. This local description removes any assumption about the orientation of the space velocity vector.

Let us examine what the notion of viscosity represents in discrete mechanics. The rotation term of the equation of motion is written

∇ d × ψ o -dt c 2 t ∇ × v .
For time constants greater than 10 -11 s the propagation of transverse waves is attenuated and the quantity dt c 2 t must be replaced by ν. However, ν does not represent a dissipation as in classical mechanics to which it is associated. The quantity ν ∇ × v represents the instantaneous rotation and the energy is really dissipated in heat only if it is not accumulated in the potential vector ψ o ; the dissipation is represented by the attenuation factor where α t ≈ 1 for fluids.

No transformation of these two compression and rotation terms is required. The first term of the second member of (10) corresponds physically to the transport of the scalar potential by the velocity field and the second one translates the dissipation of the kinetic energy by viscous effects. Indeed this energy is not accumulated during the temporal process, it is transformed into heat. It should be noted that the discrete model is full-compressible; the dissipation effects of compressional waves are translated by a factor α l different from unity and φ o would then represent the reversible part.

Assuming that viscous effects are predominant, the equation would become:

1 2 d|v| 2 dt = -v • ∇ d × (ν ∇ × v) (12)
Since the viscous term is linear, it can be modeled, as for a porous medium, by a term of the type ν/K v where K is an effective permeability.

1 2 d|v| 2 dt = - ν K |v| 2 (13)
The kinetic energy decreases therefore as E k (t) ∝ exp (-ν/K t). In the simulations the quantity-dE k /dt will be calculated directly in finite differences.

Let us consider the inviscid case ν = 0; the equation on the kinetic energy becomes:

1 2 d|v| 2 dt = -v • ∇φ (14)
In the absence of viscosity the evolution of the kinetic energy is not zero as shown by the equation ( 14) since v and ∇φ are collinear. It depends on the problem considered but, in general, the second member of this equation is the source of the divergence at infinity of the kinetic energy in direct numerical simulations without viscosity.

In the case of a decrease of the turbulence, the inertial phase can potentially present a growth of the kinetic energy which induces ineluctably the appearance of instabilities followed by a cascade of the scales of the turbulence towards the small eddies. This transfer phase ends, for viscous fluids, with the dissipation of vortices at Kolgomorov scales.

In discrete mechanics the turbulent quantities defined a posteriori are not the same as in continuum mechanics. Each quantity must be associated to a specific entity of the geometrical structure of the figure [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF], the segments for the kinetic energy, the vertices for the compression energy and the facets for the rotational energy.

E k = 1 [L] Γ * 1 2 v • v dl, E c = 1 [Ω] Ω * | φ | dv, E r = 1 [S] S * |ν ∇ × v| ds ( 15
)
where L is the length of the Γ * collection of all Γ segments of the primal geometric structure, Ω is the volume of the sum of the dual Ω * volumes associated with each vertex, and [S] is the sum of the [S * ] area of each of the S primal facets. These three quantities [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] are energies per unit mass, the kinetic energy E k , the compression energy E c and the rotation energy E r . It is important that the comparison of the turbulent quantities is performed on quantities of the same nature. The choice of the integration structure, segment, dual volume, primal surface, is adapted to the quantity concerned. For example the potential ψ, an axial vector, is defined on each facet S of the primal geometry of the figure [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF]. If the choice of volume were adopted the volume this turbulent statistic would not converge. Comparisons are more difficult in classical mechanics where kinetic energy and enstrophy are not expressed with the same units; moreover the latter is a purely kinematic quantity. These choices do not change much on the interpretation of the results of the direct numerical simulations, the calculation of the kinetic energy E r calculated on Γ (15) and that on the volume Ω follow the same evolutions in the course of time. Of course, the way the diagnostics are calculated does not change the results of simulations.

Taylor-Green vortex

The turbulence decay for the Taylor-Green vortex case must be analyzed in a more general framework where the total energy injected into a periodic cavity at the initial time is totally dissipated as heat when time tends to infinity. The initial total energy is of course composed of the kinetic energy E o k but also of the compression potential energy E o c . Indeed, the initial mechanical equilibrium must be ensured by the solution (φ o , v) and part of the energy necessary to ensure this equilibrium is associated with the compression. The energy balance can be written as follows:

E o k + E o c = E ∞ r ( 16
)
where E ∞ r is the instantaneous energy E r (T ) averaged over a period of time T tending to infinity. The equation of motion ( 8) is a law of conservation of total energy which translates the evolution of the solution of the system (φ o , v) from an initial state fixed by (E o k , E o c ) to a final state where the energy E ∞ r is in the form of heat. However, nothing prevents the transformation, at any time, of the potential energy into kinetic energy following the dynamic entanglement process.

2D Taylor-Green vortex

Before developing the results of the simulations in three dimensions of space on the TGV case it may be useful to show the consistency of the discrete physical model in two dimensions. Let us consider the domain L 2 = [-1, 1] 2 filled with a fluid of kinematic viscosity ν, the Reynolds number having then the definition Re = v o L/ν where v 0 = 1 is the maximum velocity in the domain. The initial condition imposed on the velocity is written:

   u = -v 0 cos(π x) sin(π y) v = v 0 sin(π x) cos(π y) (17)
In discrete mechanics the potential φ o (or the pressure) is not imposed, it results from the update process from ∇•v. Whatever the potential φ o imposed at the initial instant the scalar po-tential field adapts over a time constant of order of magnitude of 1/c 2 l to the value corresponding to the mechanical equilibrium.

This two-dimensional problem, notably approached by [START_REF] Sengupta | Non-linear instability analysis of the twodimensional Navier-Stokes equation: The Taylor-Green vortex problem[END_REF][START_REF] Sharma | Vorticity dynamics of the three-dimensional Taylor-Green vortex problem[END_REF], has been the subject of numerous works to exhibit turbulent behavior under certain conditions. In general the two-dimensional vortex flow is stable.

Simulations have been performed for three Reynolds numbers, Re = 100, 1000, ∞. The results are summarized by the evolution of the kinetic energy normalized by the initial value on figure (2). The evolution of the kinetic energy E k /E ko as a function of time is consistent with the fact that a 2D flow cannot be turbulent a priori if no perturbation is injected in the domain. The theory predicts that the kinetic energy decreases following an exponential law E k ∝ exp(-ν t). The case of an inviscid fluid, ν = 0, can be maintained over a long period of time with a value of the kinetic energy equal to its initial value at a precision close to that of the machine. However, for higher values of time, t > 50, the rounding errors accumulate and the solution diverges very quickly. The case ν = 0 is a solution of the problem but it is an unstable branch.

In fact, for an inviscid fluid, the imposed solution is a stationary solution, ∂v/∂t = 0, if the equation ( 18) is perfectly solved.

∇ Å |v| 2 2 ã -∇ d × Å |v| 2 2 n ã = -∇φ o (18)
It means that the inertial rotation energy is equal to the compression energy. The equation [START_REF] Ranocha | Discrete Vector Calculus and Helmholtz Hodge Decomposition for Classical Finite Difference Summation by Parts Operators[END_REF] shows that the kinetic energy E k does not vary (for ν = 0) if the term v • ∇φ is null. This is the case when the potential (pressure) field balances the inertia field. In a classical description, in the absence of variation of the pressure in a direction orthogonal to the considered plane, the pressure gradient is orthogonal to the velocity vector. This phenomenon is not related to the dimension of the space considered. The velocity field [START_REF] Ahusborde | A primal formulation for the Helmholtz decomposition[END_REF] imposed as an initial condition in a three dimensional space L 3 = [-1, 1] 3 produces exactly the same results for an inviscid fluid. If the potential (pressure) field is in initial equilibrium with the velocity field, the fields remain stationary if of course no perturbation is introduced.

3D Taylor-Green vortex

The case of the Taylor-Green vortex defined in many publications in three space dimensions is different in that the two components of the velocity (u, v) are functions of the z direction. The velocity field v = u(x, y, z) e x + v(x, y, z) e y + w(x, y, z) e z of initial 3D Taylor-Green vortex in the domain L 3 = [-π, π] 3 is of the form:

           u = -v 0 cos x sin y cos z v = v 0 sin x cos y cos z w = 0 (19)
where v 0 is a constant chosen here equal to unity.

The equilibrium pressure field deduced from the equation of motion is written:

φ o = φ 0 + v 2 0 16 (cos x + cos y) (cos z + 2) ( 20
)
where φ 0 is a constant. The simulations are performed with an initial field potential φ o = 0. As the physical model is full-compressible, waves propagate in the field with a longitudinal velocity equal to c l = 10 2 ms -1 but tests with higher velocities (10 6 ) do not change anything about the flow behavior. We simply observe that they attenuate very quickly to give the solution [START_REF] Caltagirone | On primitive formulation in fluid mechanics and fluid-structure interaction with constant piecewise properties in velocity-potentials of acceleration[END_REF] after some oscillations due to the compressible motion.

The solution of this problem has been obtained many times from various numerical methods. The development of spectral and pseudo-spectral methods notably by Orszag [START_REF] Orszag | Analytical theories of turbulence[END_REF][START_REF] Gottlieb | Numerical Analysis of Spectral Methods: Theory and Applications[END_REF] and used by many other authors [START_REF] Brachet | The Taylor-Green vortex and fully developed turbulence[END_REF] are very well adapted to simulate turbulence. The smallest structures of the turbulence are captured with a high accuracy allowing to represent the energy cascade at high Reynolds numbers. Since then, other high precision methods have been implemented to find the same results; it is for example the case of the Discontinuous Galerkin method [START_REF] Van Rees | A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical fows at high Reynolds numbers[END_REF][START_REF] Wang | High-order cfd methods: Current status and perspective[END_REF][START_REF] Debonis | Solutions of the Taylor-Green Vortex problem using high-resolution explicit finite difference methods[END_REF] which gives excellent results on the TGV case.

The reference case used as a benchmark to specify the properties of the numerical methods in direct simulation corresponds to a Reynolds number of Re = 1600. At this Reynolds number the smallest scale of viscous dissipation is indeed very small and justifies the choice of efficient methods. The results of the different authors are very close and a synthesis is available [START_REF] Van Rees | A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical fows at high Reynolds numbers[END_REF], in file "spectral-1600-512.gdiag".

For the results obtained from the Navier-Stokes model the temporal evolution of the kinetic energy is calculated on the volume:

E k = 1 [Ω] Ω ρ 2 v • v dv (21)
where [Ω] is the domain volume. The kinetic energy dissipation rate is, for an incompressible flow:

ε = 2 ν E = - dE k dt = 2 ν Ω ρ 2 ω • ω dv ( 22
)
where E is temporal evolution of the enstrophy integrated on the domain Ω.

For compressible flows other terms must be considered but as the Mach number is very low they will be ignored. The kinetic energy E k and the dissipation -dE k /dt obtained by a spectral method [START_REF] Van Rees | A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical fows at high Reynolds numbers[END_REF] with a spatial approximation of 512 3 are reproduced on the figure (3a). We observe a monotonic decay of the kinetic energy from its initial value and an always positive dissipation.

dE k dt = -2 ν E (23)
If the viscosity of the fluid is zero the global kinetic energy is conserved over time.

A direct simulation has been performed from the discrete equation of motion [START_REF] Caltagirone | Discrete Mechanics, concepts and applications[END_REF] under the same conditions as those defined by many authors for a Reynolds number of Re = 1600 and an initial condition fixed by the relations [START_REF] Lipnikov | Mimetic finite difference method[END_REF]. The results obtained are given in figure (3b). They differ significantly from those obtained by many authors with high precision numerical methods.

In spite of the different definitions of the calculated turbulent quantities a posteriori with those of discrete mechanics [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF], the calculation of the velocity and pressure are strictly corresponding and can be compared. There are noticeable differences in these first-order quantities, so there is no need to calculate the higher-order moments. The kinetic energy of the discrete mechanics model increases in the purely inertial region before starting the energy cascade towards the smallest scales. At the same time the magnitude-dE k /dt negative initially becomes positive before becoming globally constant in the course of time to then tend towards zero due to the preponderant viscous effects at this stage.

Before attempting to explain the differences observed between these two models, it is necessary to show that the results of discrete mechanics are not due to numerical artifacts.

Numerical validation

In order to confirm the growth of the kinetic energy with time during the inertial phase, simulations were performed at Reynolds numbers of Re = 100, 500, 1600 and for an inviscid fluid [START_REF] Sagaut | homogeneous turbulence dynamics[END_REF].

simulation carried out for a Reynolds number of Re = 100 shows that the flow does not have a turbulent character but the kinetic energy increases, passes through a maximum we observe an abrupt break for a time of about t ≈ 4.5 associated with the point of inflection of the curve E k (t). Thus the growth of the kinetic energy seems to be a persistent phenomenon independent of the value of the Reynolds number.

The evolutions of the kinetic energy and dissipation for Reynolds numbers of 500 and 1600 being very close, it becomes useless to continue the investigations for the Reynolds number of 1600 which requires higher spatial approximations. Indeed, while the kinetic energy evolutions are perfectly continuous, we observe on the figure (3) fluctuations on the dissipation evolutions at large time constants which are the sign of an under-resolution. The simulations are performed on a laptop computer with a processor at 2.4 GHz which limits the use of higher approximations. The objective here is to understand the significantly different results of the discrete mechanics compared to the reference model that is the Navier-Stokes equation.

The next step is to remove any doubt on the quality of the results of the simulations performed in discrete mechanics. The numerical method associated with this one is of order two in space and time; the results on many analytical solution cases (Couette flow, Poiseuille flow, ...), numerical benchamrk (steady flow past a circular cylinder, lid-driven cavity, backward facing step, twophase flows, ...), show results strictly identical to those of the Navier-Stokes equation.

Although the numerical methodology is of order two in space and time, many differences exist between the methods. The discrete equation ( 8) does not require any additional discretization, the differential operators ∇ • v, ∇ × v, ∇φ, ∇ d × ψ are used as they are geometrically defined. There is no interpolation, the physical properties, c l and ν are associated to the operators, respectively ∇ • v and ∇ × v. The sum of the operators -∇φ + ∇ d × ψ is projected on the same segment Γ. The stencil corresponding to an unknown, v associates 74 neighbors that is a linear system where each row has 75 nonzero terms; the number of unknowns is n u = 3 n e where n e is the number of segments and the number of nonzero terms in the close matrix is nnc = 75 n u . The nonlinear terms are linearized as ∇ d × (v n • v n+1 /2) where v n+1 is the solution of the linear system. The solution of this system is realized by a unpreconditioned conjugate gradient method BiCGStab2 in a few tens of iterations whatever the number of degree of freedom. This methodology gives the equation particular properties of robustness and effective accuracy.

A series of simulations are therefore undertaken to verify that the global methodology allows All quantities converge as a function of the spatial approximation, it is then possible to calculate an order of convergence. This convergence study is carried out from a Ridchardson extrapolation on E k using an error norm L 2 . The figure [START_REF] Orszag | Analytical theories of turbulence[END_REF] shows indeed a convergence at order two in space on the kinetic energy. These results corresponding to Re = 500 allow to qualify the simulations at least up to this Reynolds number to interpret them from the point of view of the described physics.

The vortex stretching phase can be observed for t ∈ [0, 3] and represented by a pressure value φ o = -0.1 on the figure [START_REF] Gottlieb | Numerical Analysis of Spectral Methods: Theory and Applications[END_REF] for a Reynolds number of Re = 1600. Given the closed geometry of the TGV case, the elongation of the vortices is not possible but the reduction of the isobaric section with time is very marked. Consequently the rotation velocity, fixed by the scalar φ i , increases significantly which translates into an important increase of the kinetic energy E k .

The evolution of vortex structures is much better discernible on Bernoulli φ o B scalar potential fields. The figure [START_REF] Caltagirone | Discrete Mechanics, concepts and applications[END_REF] shows some snapshots of isovals of the Bernoulli scalar potential φ o B superimposed on the inertial potential field φ i = |v| 2 /2 for times T = 0, 1, 5, 20 for a Reynolds number of Re = 1600. The visual results of the velocity or scalar potential fields are quite comparable to those of the literature on the same subject [START_REF] Van Rees | A comparison of vortex and pseudo-spectral methods for the simulation of periodic vortical fows at high Reynolds numbers[END_REF][START_REF] Diosady | Case 3.3: Taylor-green vortex evolution, Case summary for 3rd International Workshop on Higher-Order CFD Methods[END_REF].

The initial condition on φ o is set to 0 because this quantity must be subject to the velocity but after a very short time the solution on φ corresponds to the pressure value given by [START_REF] Caltagirone | On primitive formulation in fluid mechanics and fluid-structure interaction with constant piecewise properties in velocity-potentials of acceleration[END_REF]. The pressure field for t = t b = 1 reveals the beginning of the vortex bursting process preceded by a vortex-stretching phase where the local velocity increases, the vortex radius reduces and the pressure increases. For a time t = 3 we observe the vortex bursting phase where ejections are emitted from the main vortices. The secondary flows are organized into intermingled vortices (t = 5) of intermediate sizes. For larger times all structures seem to emit smaller vortices. This scenario is the one described to explain the energy cascade in the transfer phase is similar to what is described in the literature. Despite the very different kinetic energy evolutions from those observed with the Navier-Stokes model in the inertial zone, the degeneracy of the vortices to smaller scales seems very comparable with previous results.

The analysis of the results of the evolution of the kinetic energy E k , of -dE k /dt and of the potential φ o as a function of time allows us to find the behavior inscribed within the discrete kinetic energy theorem for the inertial phase (10):

1 2 d|v| 2 dt = -v • ∇φ o (24)
Indeed the increasing evolution of -dE k /dt is in t 2 and the decreasing one of φ o is also in t 2 . The calculation of the time evolution of v shows that the velocity is constant on average. Thus the kinetic energy varies as E k ∝ t 3 .

The figure [START_REF] Caltagirone | An alternative to the concept of continuous medium[END_REF] traces schematically the variations of the kinetic energy for the inertial phase, the transfer phase and the dissipation phase. From t = 0 to t = t s the inertial phase is accompanied by a dissipation even if it is relatively weak as shown in figure [START_REF] Sagaut | homogeneous turbulence dynamics[END_REF]. The first burst appears for a time t b < t s but the transfer phase really starts at time t s after smaller vortices interact. The evolution of -dE k /dt shows that in the energy cascade this quantity fluctuates but remains almost constant on average hence the approximation in t. In the dissipation zone, the quantitative examination of the evolution of the kinetic energy E k (t) in log-log coordinates shows an exponential law E k ∝ exp(-ν t).

From the energetic point of view the growth of E k is closely associated with the decrease of the compression energy E c ; in this phase the angular momentum is indeed conserved. Recall that there is no injection of kinetic energy, the domain is periodic in all directions and no source term is imposed. This energy transfer is essentially inscribed in the entanglement between inertia, compression and rotation effects. In the absence of significant viscous shear the compressional energy is converted to rotational energy, which is what is meant by the vector potential ψ o = ν ∇ × v which also carries the block rotation of the vortices. Energy is truly degraded into heat only if the spatial scales are small enough to be dissipated by velocity gradients.

The spectral analysis of the results obtained is for the moment put aside because it requires simulations at higher Reynolds numbers to significantly recover the undeniable results on turbulence, in particular the -5/3 slope of the Kolgomorov theory. The attention is rather focused on the physical analysis of the behavior of the Navier-Stokes equations and the discrete mechanics.

Discrete turbulence analysis

It should be recalled that the discrete model ( 8) has allowed to strictly recover the results of the Navier-Stokes equation for fluid flows and those of the Navier-Lamé equation for solids and their coupling for Fluid-Structure Interaction. Turbulence is addressed for the first time with this physical model. The conclusion of the previous section concludes that the attribution of significant differences on the TGV case simulations is only due to the physical model, discrete mechanics and Navier-Stokes. The differences are now analyzed in depth from the properties of the Navier-Stokes equation and the discrete equation of motion.

Entanglement of compression and rotation terms

The strong coupling between all the terms of the discrete equation of motion is essential for its physical representativity. Indeed, the application of the divergence or primal curl operators on the equation of motion, for example the Navier-Stokes equation, modifies its meaning. The divergence operator is used in the time-splitting methods to project the equation on a zero divergence field and the application of the primal curl operator leads to the elimination of the pressure term. These procedures break the entanglement between the compression and rotation terms. In most cases this is of little importance, especially in the steady state, but the treatment of turbulence becomes problematic when operators are used to transform the equation of motion.

incompressibility constraint ∇ • v = 0 is then imposed without any possibility of transfer by compressional waves.

For example the application of the primal curl to the incompressible Navier-Stokes equation allows to eliminate the pressure gradient -∇p to make appear the vorticity form:

     dω dt -ω • ∇v = ν ∇ 2 ω ∇ 2 Ψ = -ω (25)
where ω = ∇ × v is the vortex vector and Ψ is the potential velocity vector. Many methods use this formulation which has some advantages including the fact that it ensures the incompressibility constraint ∇ • v = 0. However, the formulation (25) inhibits any transformation between the compression and rotation energies. Incompressibility is a very restrictive assumption insofar as it propagates any pressure perturbation instantaneously to infinity. This quantity becomes a simple Lagrangian to ensure the mass conservation constraint under the incompressibility assumption.

The discrete primitive formulation represented by the equation ( 26) has a remarkable property due to the Helmholtz-Hodge decomposition.

∂v ∂t = -∇ φ o + φ i -c 2 l dt ∇ • v + ∇ d × ψ o + ψ i -c 2 t dt ∇ × v (26)
The two terms of the right-hand side are locally orthogonal [START_REF] Caltagirone | On Helmholtz-Hodge decomposition of inertia on a discrete local frame of reference[END_REF] but they are carried by the same Γ support. The compression energy represented by the first term cannot exchange directly with the last term which translates the rotation. This transfer is only possible if the motion is unsteady; the energy is thus redistributed towards one or the other term of the right-hand side. This is the remarkable idea of J.C. Maxwell's of combining magnetism and electrodynamics into electromagnetism. These two effects are dynamically intertwined.

The reasons for the breaking can be multiple, the splitting of the operators can be one of them. The Navier-Stokes equation translates the conservation of momentum and the law of conservation of mass is attached to it, the Navier-Stokes equation does not ensure as such the conservation of mass. In an incompressible formulation, the pressure is deduced from this last law by assuming a link, through a state law, between the density and the pressure. Numerical methods such as the projection method [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF] add an additional splitting.

It is quite different for the equation ( 26), even for motions considered incompressible. Indeed the term dφ = -dt c 2 l ∇ • v never becomes zero, when the celerity c l increases it is the divergence that decreases. The quantity dφ is an energy per unit of mass transferred to the fluid by a compression process during the time lapse dt. The dynamic entanglement takes all its meaning in turbulence by definition unsteady and where the frequencies are high.

The law of discrete motion (8) expresses primarily the conservation of total energy per unit mass. The left term represents the variation of the total mechanical energy between a and b and the two right terms the compression and rotation energies; each of the last two terms is a Lagrangian with a potential energy fixed by the potentials φ o and ψ o , and a kinetic energy carried by the increments dφ = dt c 2 l ∇ • v and dψ = ν ∇ × v. More precisely the equation of motion can be put in an integral form carried by the segment Γ:

a ∂v ∂t • t dl = - b a ∇ Å φ o -dt c 2 l ∇ • v + |v| 2 2 ã • t dl + b a ∇ d × Å -ν ∇ × v + |v| 2 2 n ã • t dl (27)
where ∂v/∂t represents, for an averaged flow, the velocity fluctuation, ∇ dt c 2 l ∇ • v the acceleration which represents the compression fluctuation and finally ∇ d × (ν ∇ × v) the acceleration associated with the rotation fluctuation; this operator filters the local block rotation eliminated by the dual curl.

In the Helmholtz-Hodge decomposition

∂v ∂t = -∇ (φ o + φ i ) + ∇ d × (ψ o + ψ i ) (28)
the total energy represented by the integral on a segment of the first member cannot be dissociated from one of the components, it is the total energy that is conserved and this is the very object of the law of dynamics. At each instant the compression energy and the shear energy are exchanged only through the term in time ∂v/∂t because these two terms are orthogonal and exchange nothing directly. This concept is called dynamic entanglement; it comes from the remarkable idea of J.C. Maxwell who founded electromagnetism by combining the laws of electrodynamics with those of magnetism. An alternating current in an electrical conductor produces an induced electric current of the same frequency in a circuit surrounding the primary conductor, if the current is direct the coupling does not exist; it is this observation that led to the success of alternating current for the distribution of electricity.

In mechanics, this concept translates into the circulation of a fluid in a circuit composed of Γ segments such as those forming the primary surface of the figure [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF]. This circulation allows to calculate the primal curl to define the potential vector ψ = ν ∇ × v. In its turn the circulation of this dual vector reprojects a flow on the segment Γ. The two accelerations -∇φ and ∇ d × ψ circulate on the same segment without interacting, the two components are indeed orthogonal. On the other hand, one of these components modifies the velocity in time, which has the effect of affecting the other one a posteriori.

In turbulence this mechanism of transferring compression and rotation energies is essential. In classical mechanics the role of compression is underestimated if not neglected. But even for flows considered incompressible, it is not the fact that ∇ • v is close to zero that is important, what really counts is the quantity dt c 2 l ∇ • v which remains defined when ∇ • v = 0 because the celerity c l tends to infinity.

A local point of view

Continuum mechanics is part of the larger concept of continuous medium where any quantity is first defined on a bounded domain Ω limited by a surface Σ before being localized at a point from the divergence theorem. Thus the momentum and mass balances are first established on the Ω domain and then the different quantities are assigned at a point. Let's consider the momentum equation associated to the conservation of mass:

         ρ dv dt = ∇ • σ + f dρ dt + ρ ∇ • v = 0 ( 29
)
where σ is the Cauchy stress tensor and f a source term per unit mass. Boundary conditions complete the local formulation. kinetic energy theorem obtained by multiplying scalar by v and integrating over the volume can be interpreted as a kinetic energy balance equation:

d dt Ω ρ |v| 2 2 dv = Σ σ v • n ds - Ω σ : D dv + Ω v • f dv ( 30
)
where D is the strain rate tensor. For a Newtonian fluid with constant physical properties such that σ = -pI+τ where τ is the viscous stress tensor and leaving aside the external forces, the kinetic energy theorem becomes:

d dt Ω ρ |v| 2 2 dv = - Σ p v • n dv + Ω p ∇ • v dv + Σ v τ • n ds - Ω τ : ∇v dv (31)
When the velocity is zero or v • n = 0 at the edge the first integral of the second member disappears, the second integral also disappears if the flow is incompressible as well as the third if the normal stresses are zero at the edge; this is the case for the Taylor-Green vortex example. The kinetic energy depends only on the viscous dissipation term represented by the last integral. For an inviscid fluid the kinetic energy is conserved during the motion. The elimination of the pressure from the kinetic energy equation [START_REF] Ranocha | Entropy Conserving and Kinetic Energy Preserving Numerical Methods for the Euler Equations Using Summation-by-Parts Operators[END_REF] used in the same way for the solution of the Navier-Stokes equation by some variational formulations leads to weak solutions. See [START_REF] Ranocha | Entropy Conserving and Kinetic Energy Preserving Numerical Methods for the Euler Equations Using Summation-by-Parts Operators[END_REF] for a relevant discussion on kinetic energy preservation for the Euler equations.

It is quite different in discrete mechanics, this formulation corresponds to a strong solution of the equation of motion. The solution is sought on each of the Γ segments of the primal structure. The conservation is expressed as all other accelerations on this segment. The continuous medium frame is completely abandoned and replaced by a local discrete frame. The classical Galilean frame of reference is replaced by the local frame of reference of the figure [START_REF] Richardson | Weather Prediction by Numerical Process[END_REF]. As v is the velocity component on Γ and ∇φ is the restriction of the potential gradient (pressure), these two vectors are collinear.

dE k dt = -v • ∇φ -v • ∇ d × (ν ∇ × v) (32)
The first term of the second member -v •∇φ is not zero a priori and its sign is indeterminate. The last term represents the viscosity effects whose contribution is definitely converted into heat, it is the viscous dissipation. It is possible to model this term from Darcy's model by introducing an effective permeability K to transform it into -(ν/K) |v| 2 . In any case, since ν and K are positive quantities, this term tends to reduce the kinetic energy with time following an exponential evolution.

The turbulent statistic -dE k /dt called dissipation in continuum mechanics is a by-product of the numerical solution. It does not depend fundamentally on the way it is computed, on the volume in continuum mechanics or on a segment in discrete mechanics. However, it reveals that the kinetic energy calculated during the resolution has a significantly different evolution by the two approaches.

The local point of view adopted in discrete mechanics makes it possible to ignore long distance influences, in particular those of the boundary conditions. In flows with very small time constants dt like those of turbulence, the interactions are limited in space by the discrete horizon dh = c l dt, a fluctuation of a quantity at a point is not instantaneously felt by the boundaries of the domain. The interactions from one segment to another are actualized from cause to effect through the scalar potential common to two segments. Of course the modeling of a viscous flow with a larger time constant forces to consider the transformation of the transverse propagation term

-dt c 2 t ∇ × v into a diffusion term -ν ∇ × v.
This assumption is justified by the very rapid decay of transverse waves into heat by viscous dissipation.

The role of pressure is largely neglected in some models of turbulence, which introduce a turbulent pressure term that is associated with viscosity. The notions of pressure and viscosity are strictly disjoint in discrete mechanics. The pressure (scalar potential) is also largely underestimated in direct simulations performed in incompressible or compressible formulations using constitutive laws. Turbulence is characterized by a very small scale locality and violent phenomena. The Taylor-Green vortex phenomenology is composed of a succession of vortex-stretching, collapses and vortex bursting. The figure (10) shows the time evolution of the mean pressure φ o and its maxima in the domain.

We observe that the minimum pressure increases at first and then abruptly changes its behavior for a time t b ≈ 3; this time corresponds to the end of the vortex-stretching phase of the four vortices of axis z. We observe that the local pressure variations are of order of magnitude |v| 2 , values much larger than the mean pressure. The next phase is a rapid decrease of the pressure which stops for a time t ≈ 5 followed by multiple relaxations and compression. The evolutions of the pressure correspond to faster and faster frequencies denoting a reduction in the size of the vortices generated by the vortex bursting. For this Reynolds number of Re = 500 the effects of the viscosity become preponderant for a time t v ≈ 15. In fact the transfer phase does not correspond to a regular decrease of the kinetic energy; indeed the structure of the flow is marked for a long time by the initial condition. On the evolution of -dE k /dt we can observe the rise of this quantity, which corresponds to local collapse phenomena followed by vortex-bursting. Obviously the evolution of the kinetic energy is more regular; it allows to quantify the law of decrease of E k as a function of time. This transfer phase which ends approximately at a time t v is correlated to the change of behavior on all quantities, E k , -dE k /dt, E c , E r . Compression and rotation effects are always present and of the same order of magnitude for a fully developed turbulent flow. Each of these two effects is preponderant respectively in the inertial zone and in the viscous dissipation zone. Their entanglement is of course prevalent in the transfer zone from large vortices to small vortex structures, hence the importance of the mathematical formulation adopted to perform direct simulations.

On the increase of kinetic energy

The quantity -dE k /dt is classically called kinetic energy dissipation rate. In fact, this quantity is not only linked to the dissipation. The equation (10) shows that this quantity can vary independently of the viscosity even when ν = 0. In a viscous or non-viscous flow the kinetic energy is not conserved. The equation ( 8) reflects the preservation of total energy but not of kinetic or potential energy.

No principle justifies a priori that the kinetic energy must decrease monotonically for an inviscid fluid. Let us consider a rotating tube of radius R and length L whose velocity field is given by the velocity components [START_REF] Lipnikov | Mimetic finite difference method[END_REF]. The dynamics of this tube is fixed by the conservation of angular momentum expressed in kg m 2 s -1 whose transposition in discrete mechanics would be expressed in m 2 s -1 ; it is the product of the velocity v by the distance between the material and the rotation axis. This notion of angular momentum L o = OM ∧ v is too restrictive because it is associated with a given velocity. The potential vector ψ = ν ∇ × v generalizes this concept by introducing the quantity ν which transforms it into energy per unit mass. Thus when the rotation tube lengthens the potential vector is conserved as its volume and the radius decreases. The velocity v increases as v ≈ 1/R and the kinetic energy varies as E k ≈ 1/R 2 . The mechanical equilibrium between the effects of compression and rotation can only be assured if the pressure decreases according to a law φ o ≈ t 2 . This instability continues, the rotation velocity increases, the radius R tends to zero and the pressure suddenly increases significantly. The current tube collapses locally and literally explodes producing a vortex-bursting. From this moment on, smaller vortexes with an orthogonal axis to the main vortex form in a complex tangle. This phenomenon was described in several publications in particular by Spalart [Spa98] or [Moe05]. More recently the simulation of a single vortex by van Rees in 2020 shows well the mechanism of formation of this phenomenon [START_REF] Van Rees | Vortex bursting[END_REF].

It seems that all simulations on the TGV case show a monotonic decrease of the kinetic including in the inertial zone. The figure (4) reveals a reduced effect of the viscosity in this zone but does not modify the growth character of E k with time. This is in discrete mechanics a robust mechanism. The attribution of this essential difference on the evolution of E k in the inertial zone is thus due to the choice of the physical model, Navier-Stokes or discrete mechanics. It would be risky to draw conclusions about the representativeness of the discrete model for turbulence. It must be noted that both models show exactly the same solutions for classical flows [START_REF] Caltagirone | On primitive formulation in fluid mechanics and fluid-structure interaction with constant piecewise properties in velocity-potentials of acceleration[END_REF][START_REF] Caltagirone | On Helmholtz-Hodge decomposition of inertia on a discrete local frame of reference[END_REF][START_REF] Caltagirone | Application of discrete mechanics model to jump conditions in two-phase flows[END_REF].

However, it is possible to attempt an explanation based on the two equations because they are significantly different. The discrete equation ( 8) intrinsically ensures the dynamic entanglement between the viscous and rotational effects. This is not the case for the incompressible Navier-Stokes equation because the diffusion of momentum is written as ν

∇ 2 v = ν ∇(∇•v)-ν ∇×∇×v;
as the incompressible model must satisfy the constraint ∇ • v = 0 the first term is taken equal to zero without precaution. One must also take into account the term existing in the compressible formulation λ∇(∇•v) for which the value of the compression viscosity λ for a fluid is not defined. The Navier-Stokes equation does not intrinsically satisfy the conservation of mass, it must be added to a specific equation dρ/dt + ρ ∇ • v = 0 or directly ∇ • v = 0. This decoupling inhibits entanglement phenomenon. The treatment in compressible formulation of the Navier-Stokes equation does not change the result of the monotonic decay of the kinetic energy, it introduces another form of decoupling by integrating an additional state law. Moreover, the majority of simulations performed from the Navier-Stokes equation use a splitting where the pressure is calculated by a separate Poisson equation.

Contrary to classical mechanics where the equations of motion deal with the balance of momentum and momentum, the equation ( 8) contains both notions in a form that translates the conservation of energy per unit mass. The acceleration γ carries both the conservation of motion in one direction and that represented classically by the angular momentum of classical mechanics.

Curvature in space and time of inertial potential

Another explanation can be given for the differences between the continuous and discrete formulations, it concerns the formulation of inertia. In classical mechanics this is written indifferently v • ∇v or ∇ • (v ⊗ v)v ∇ • v or ∇(|v| 2 /2)v × ∇ × v. The application of the differential operators divergence and curl bring out terms whose meaning is uncertain [START_REF] Caltagirone | On primitive formulation in fluid mechanics and fluid-structure interaction with constant piecewise properties in velocity-potentials of acceleration[END_REF]. In particular, the divergence of the Lamb vector L -v × ∇ × v produces a term which is the second invariant of the tensor ∇v which appears as a compatibility condition for solving an incompressible flow. The inertia is written in mechanics in the generic form of a Helmholtz-Hodge decomposition; its divergence and its curl eliminate from the start the dual rotational and the gradient respectively.

In a certain number of simple flows (Couette, Poiseuille, ...) the inertial effects are apparently absent. Indeed, for example for the Poiseuille flow the term v • ∇v is identically zero. But the term ∇(|v| 2 /2) is not and it is therefore necessary that the second component of the inertia is equal and opposite to the first. This is the case in continuum mechanics where the Lamb vector is written as L = -v × ∇ × v. As the classical formulation is linked to the existence of a global Galilean reference frame, the equality of the absolute values of the two terms imposes that the Lamb vector is also a gradient. This is not the case in discrete mechanics where the gradient of a function can be equal or opposite to the dual curl of a different function. The inertia term is therefore also zero by compensation in discrete mechanics. The material derivative is then written in the form:

dv dt = ∂v ∂t + ∇ Å 1 |v| 2 ã -∇ d × Å 1 2 |v| 2 n ã (33)
The properties ∇ × ∇φ = 0 and ∇ • ψ = 0, strictly verified whatever the primal and dual structures, allow to filter the equation of motion if necessary in the framework of a decoupling cancelling of course the dynamical entanglement.

sum of the two components of the inertia noted κ i = ∇φ i -∇ d × ψ i can be interpreted as the average curvature of the inertial potential φ i = |v| 2 /2 [START_REF] Caltagirone | On Helmholtz-Hodge decomposition of inertia on a discrete local frame of reference[END_REF]. The analogy between the calculation of the mean curvature as the sum of the two local curvatures of a surface in two orthogonal directions tangent to it seems indeed relevant. In discrete mechanics all physical effects are associated with the Helmholtz-Hodge decomposition in the form of two components, one curl-free and the other divergence-free. The curvature in time of the kinetic energy E k visible on the figure (4) for ν = 0 adds a complementary notion. The increase of E k with time, approximately in ∝ t 3 , takes place with a strong negative curvature. Furthermore, the gravitational potential per unit mass is written φ g = -G M/r where G is the gravitational constant, M is the mass of the body and r is the distance of the point considered from the center of mass. Its curvature is equal to κ g = -∇φ g + ∇ d × φ g n. By analogy, if |v| 2 /2 is the inertial potential, κ i is its curvature in space and dE k /dt its representation in time.

Conclusions

Direct numerical simulations based on the Navier-Stokes equation using high precision methods implemented by many authors provide the same results. All these results acquired on the Taylor-Green vortex case are not contestable. On the other hand, the discrete equation of motion allows to find the results of classical mechanics on many analytical and numerical solutions. Moreover, the solutions of the discrete equation of motion converge at order two in space and time including the presented case. For the first time the solutions of the discrete equation of motion differ from those of Navier-Stokes on the turbulent TGV flow. The difference can therefore be attributed solely to the physical model.

Discrete mechanics is an alternative to the notion of continuous media; it is not only a proposal to modify the equations of classical mechanics. It addresses a physical model based on a discrete geometrical concept; it abandons the global Galilean frame of reference in favor of a local frame of reference where the interactions are of cause and effect but also the derivation at a point, the integration, etc. The discrete kinetic energy theorem deduced from the equation of motion brings to light a term v • ∇φ related to the advection of the potential by the velocity component; these two vectors carried by the same segment are collinear. It is this term that translates the vortex-stretching at the origin of the turbulence triggering by a vortex bursting in the case of the Taylor-Green vortex. Despite the significant difference in behavior observed in the inertial phase of the kinetic energy evolution, the energy cascade that follows, illustrated by the vortex splitting, is very close to that simulated by the Navier-Stokes equation.

This new physical model has different properties from the standard model represented by the Navier-Stokes equation. It is not in contradiction with the general principles of mechanics but presents an original discrete vision. This study, like the previous ones, is one of the steps to confirm or deny its legitimacy.

Figure 1 .

 1 Figure 1. A primitive planar facet S defined by a contour Γ * is oriented according to the normal unit n such that n • t = 0; the dual surface ∆ connecting the centroids of the cells defined by the dual contour δ is also flat.

Figure 2 .

 2 Figure 2. 2D Taylor-Green vortex for a spatial approximation 128 2 and for two values of Reynolds numbers Re = 100, 1000 and for an inviscid fluid.

Figure 3 .

 3 Figure 3. Taylor-Green vortex at Re = 1600; in top (a) the spectral Re-1600-512.gdiag simulation provides the evolution of dimensionless values of E k , ε = -dE k /dt for the Navier-Stokes model [29] and in bottom (b) the same results for the discrete mechanics model.

Figure 4 .

 4 Figure 4. Re = 500, 1600 and for an invicid fluid

Figure 5 .

 5 Figure 5. Re = 500 (a) kinetic energy E k , (b) -dE k /dt, (c) energy of compression E c and (d) energy of rotation E r .

Figure 6 .

 6 Figure 6. Convergence in space on the kinetic energy in L 2 norm for Re = 500.

Figure 7 .

 7 Figure 7. Vortex-stretching phase for Taylor-Green vortex at Re = 1600; snapshots of value of potential φ o = -0.1 colored by the inertial potential φ i = |v| 2 /2 for t = 0 and t = 0.8.

Figure 8 .

 8 Figure 8. Taylor-Green vortex for Re = 1600 in 196 3 ; snapshots of iso-values of the φ o B potential colored by the inertial potential φ i = |v| 2 /2 for t = 0, t = 1, t = 5 and t = 20.

9 .

 9 Scheme of the evolution of the kinetic energy in the inertial phases, the energy cascade and the dissipation zone.

  b

Figure 10 .

 10 Figure 10. Evolution of the maximum and minimum pressure in time at Re = 500.

These two components are assigned to the same segment Γ as the acceleration itself, which is its sum. It is essential to note that these two components have no direct action on each other, they simply circulate on the same trajectory. The Helmholtz-Hodge decomposition is sometimes written with a third harmonic component which is both curl-free and divergence-free; in the present case, this last term is deleted because it corresponds to the uniform movements of translation and of rotation separated from the outset by the discrete differential operators. The energy conservation balance between a and b is then written in the form of a discrete integral: Γ