Jean-Paul Caltagirone
email: calta@ipb.fr

A primitive formulation (v, ψ) for incompressible flows based on velocity and vector potential of the acceleration

Keywords: Navier-Stokes equation, Discrete Mechanics, Helmholtz-Hodge Decomposition, Conservation of Acceleration, Noether Theorem

An original primitive formulation of the equation of motion for viscous flows, based on two variables, the velocity and the vector potential of the acceleration, ensures the strict incompressibility of the flow. This equation, which reflects the conservation of the total energy, states that the intrinsic acceleration of the material medium is equal to the sum of the two terms of a Helmholtz-Hodge decomposition, one divergence-free and the second curl-free; the velocity field is itself decomposed into a solenoidal and an irrotational part. The incompressibility approximation of a flow leads to retain only the solenoidal components of all the terms of the equation of motion, viscous, inertial and vector potential of the acceleration; thus, all the terms which are expressed in the form of a dual curl are free of divergence. If, at a given moment, the velocity field is solenoidal, it will remain so both on a trajectory and in the whole domain. This decomposition allows us to state dual Bernoulli's theorem which applies to rotational flows as the equivalent of Bernoulli's theorem associated to irrotational flows.

Introduction

There are many formulations of the equations of fluid mechanics and each has advantages and disadvantages that are described in many textbooks in this area of physics. They apply to a variety of equations, but the most important of these is the Navier-Stokes equation, which is defined as a partial differential equation with a mass conservation equation and possibly an energy conservation equation. The system of equations is closed by constitutive laws.

For simplicity, two families of formulations are used, the formulations in primitive variables and those where a derivation operator is applied to the Navier-Stokes equation. The first one is symbolized by the couple (p, V), the pressure and the velocity. It is a natural formulation that has many qualities and is legitimately used to simulate all types of flows. A large number of computational methodologies are based on this form, the projection methods [START_REF] Chorin | A numerical method for solving incompressible viscous flow problems[END_REF][START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF], the artificial compressibility method [START_REF] Peyret | Computational Methods for Fluid Flow[END_REF] or its implicit form, the augmented Lagrangian method [START_REF] Fortin | Méthodes de lagrangien augmenté, Application à la résolution numérique de problèmes aux limites[END_REF], to mention only these techniques.

The second family of formulations of the Navier-Stokes equations is based on the application of a differential operator to conserve or eliminate the pressure. The divergence of the equation leads to a pressure equation in the form of a Poisson equation with a relatively complex righthand side. The curl operator is applied if the pressure has to be eliminated; this is the case for the rotational-stramfunction (ω, ψ) formulation in two space dimensions or (Ω, Ψ) in rotationalvector potential in three dimensions where velocity equal to V = ∇ × Ψ. This second family of methods has serious drawbacks, in particular when the physical properties are not constant. Moreover, it leads, in dimension three, to increase the number of equations to solve. However, it has the advantage of strictly safeguarding the incompressibility in the case where it is necessary.

The problem posed here is to propose a formulation which preserves the interests of a formulation in primitive variables and those of the second family where the constraint ∇ • V = 0 is strictly maintained naturally. In the framework of classical mechanics this question does not seem to arise because it is the pressure p, representative of the effects of compressibility, which maintains the mechanical equilibrium; it is the scalar potential of the intrinsic acceleration γ ≡ dV/dt or rather the ratio φ = p/ρ where ρ is the density of the medium. The equivalent of the scalar potential is the acceleration vector potential, ψ; it represents for viscous terms what pressure is for compression effects. In fact, for incompressible flows, the pressure is only a Lagrangian with no real physical existence, it is only used to reach the mechanical equilibrium and to maintain an acceptable level of divergence during a simulation. The discrete mechanics equation, whose basics are briefly recalled in the following sections, states that the proper acceleration of a material medium is equal to the sum of the gradient of the scalar potential φ and the dual curl of the vector potential ψ. From then on, the solution of the problem consists in replacing the scalar potential -∇φ of acceleration by its dual curl ∇ ⊗ ψ, these two quantities of different variables being equal to a constant. Thus, the identity ∇ • (∇ ⊗ ψ) = 0 ensures strict incompressibility.

All formulations in primitive variables containing the pressure used to ensure the incompressibility of the flow are confronted with this paradoxical difficulty, one must know a non-zero divergence at any time of the solution to obtain ∇ • V = 0. From a practical point of view, obtaining the pressure in the codes for solving the Navier-Stokes equation constitutes a very important part of the total solution time.

The starting point of the presentation devoted to incompressible flows is Newton's second dynamic law in discrete form which states that the acceleration is equal to the sum of the accelerations corresponding to the compression and rotation accelerations γ = -∇φ+∇⊗ψ. The attribution of a particular status by Galileo to mass (density) is not a necessity, the Equivalence Principle is only linked to accelerations, that of the proper inertia of the medium and that of gravitation. Newton's second law is thus transformed into a balance of accelerations.

General framework

Primal and dual structures

An important distinction must first be made about the notion of velocity; V is the velocity vector of R 3 of continuum mechanics whereas v is the scalar assigned to an oriented segment of discrete media mechanics. Similarly, one should not confuse V = ∇ × Ψ where Ψ is the vector potential of velocity and that of acceleration ψ o = ν ∇ × v where ν is the kinematic viscosity.

The notion of discrete medium allows us to introduce the basic elements of differential geometry which are used to define the interactions between the different quantities of physics. The physical domain is composed of two structures called primal and dual which can be interchanged a priori, they are represented schematically on Figure 1. The primal geometrical structure of Figure 1(a) is composed of rectilinear segments Γ oriented by a unit vector t; the extremities (a, b) of each segment define its length dh named discrete horizon on which is fixed the scalar potential of the acceleration φ. The dual structure of Figure 1(b) is formed by the closed contour ∆ which carries the vector potential of the acceleration positively oriented by the unit vector n. The primal contour Γ can be open or closed and when the family Γ * of segments Γ forms a closed circuit, the two structures are nested. Thus, the primal contour can be open or closed but the dual contour ∆ is necessarily closed. The Stokes theorem allows then to transform a vector circulation on a contour into a curl on the surface bounded by Γ; a second application of this theorem on ∆ allows to project the result on the initial contour. This discrete approach can be easily transformed into a local approach by considering a homothetic reduction, dh → 0, of the structures where the angles are conserved. Contrary to the notion of continuous medium, the result of this reduction is not a point. The local reference frame created by the assembly of the two structures does not require the creation of the global reference frame of continuum mechanics, the classical reference frame R 3 (x, y, z) does not exist. The interactions between the different quantities are realized according to the principle of causality adopted in the present context.

The discrete operators for transforming a vector into its dual or inverse are (i) the divergence of velocity, ∇ • v which assigns the result to a vertex of the primal structure, (ii) the primal curl ∇ × v which fixes the result on the barycenter of the primal facets, (iii) the gradient of the scalar is represented by ∇φ, a quantity that is uniquely defined on the segment Γ oriented by the unit vector t and finally (iv) the dual curl ∇ ⊗ ψ that reprojects the result of the circulation on the dual contour ∆ on the segment Γ. This dual curl operator does not exist in the context of the continuous medium, it is replaced by the nature of the vector to which it applies, a polar vector or an axial vector still called pseudo-vector. Several notations exist in the field of differential geometry in index or exponent of the operator ∇ but, given the importance of this operator, it is appropriate to give it a symbol which completes in a coherent way the first three definitions, it is thus the symbol ∇⊗ which is chosen. This means that the curl is computed directly by the circulation of a vector on the closed contour ∆ by the Stokes theorem. Classically, in continuum mechanics, it represents the tensor product, but in discrete mechanics, tensors do not exist and the vectors themselves are scalars on oriented segments. From a physical point of view, the entanglement of the two geometrical structures of the Figure 1 The properties of this geometric description are remarkable, in particular the classical identities of the continuous medium,

∇ • (∇ ⊗ ψ) = 0, ∇ × (∇φ) = 0, (1)
are verified exactly whatever the choice of the polygonal or polyhedral planar surfaces which tessellate the physical domain. In a continuous medium, these relations are ∇ • (∇ × V) = 0 and ∇ × (∇f ) = 0 where V and f are differentiable vector and scalar functions respectively. If these relations are mathematically correct, they have a physical representation flaw. In discrete mechanics, the divergence of a curl applies only to the vector potential ψ and the curl of gradient applies only to the scalar potential φ. If v is the vector carried by the segment Γ we get the following expressions:

∇ • (∇ ⊗ (∇ × v)) = 0; ∇ × (∇ (∇ • v)) = 0 (2)
Indeed, the scalar potential is characterized by the divergence of the velocity and the vector potential is defined from the curl of the velocity. The scalar and vector potentials of the acceleration vector are at the same level and are expressed in the same units. In continuum mechanics and electromagnetism, there is still a confusion between the potential acceleration vector and the potential velocity vector. To obtain a velocity field at zero divergence, it seems natural to pose V = ∇ × Ψ such that ∇ • ∇ × Ψ = 0 where Ψ is the velocity vector potential. Similarly, in electromagnetism, the scalar potential of the acceleration is, e, the electric potential but the vector potential A is not that of the acceleration but that of the induction field B. The induction field is conservative if we fix B = ∇ × A where A is the vector potential of the induction. This additional derivation to obtain a field with zero divergence leads to harmful and useless effects both in fluid mechanics and in electromagnetism.

In fluid mechanics, the condition ∇ • V = 0 is not an absolute necessity, it is only an approximation. The vector potential of the velocity ∇ × Ψ exists only for a strictly solenoidal field. This is not the case for the vector potential of the acceleration ψ which exists under all circumstances for any flow. The potential ψ is constructed from ∇ × v which has the effect of making the divergence-free component of the velocity of v disappear and ψ is therefore a conservative field. Similarly, the scalar potential φ is defined by ∇ • v which removes the irrotational component of v. The intrinsic acceleration then takes a structure of the form γ

= -∇ (∇ • v) + ∇ ⊗ (∇ × v).
According to the chosen approximation, some terms are naturally removed.

The law of motion is an evolution equation in time that uniquely addresses a time interval dt between t o which defines mechanical equilibrium at the previous instant and a current time t = t o + dt where mechanical equilibrium is also required. The time interval is chosen so that the discrete horizon dh is travelled at a velocity equal to or less than celerity of medium such that dt ≤ h/c. At time t o , the velocity is equal to v o and the potentials at φ o and ψ o are called retarded potentials [START_REF] Liénard | Champ électrique et magnétique produit par une charge électrique concentrée en un point et animée d'un mouvement quelconque[END_REF] in reference to those in electromagnetism. The velocity v o is a symbolic notion because the integration in time can be realized by a scheme of any order.

Equation of motion

The law of motion in terms of accelerations developed in the last few years is presented as an alternative to the Navier-Stokes equation; it allows to simulate flows of any nature by strictly finding the classical solutions of fluid mechanics [START_REF] Caltagirone | Application of discrete mechanics model to jump conditions in two-phase flows[END_REF]. Several important differences distinguish the proposed formulation from the Navier-Stokes equation. First of all, the concept of continuous medium is abandoned, as well as the one-point derivation, the integration and, in general, the mathematical analysis. The concept of a global frame of reference is replaced by the existence of local frames of reference where interactions are realized from cause to effect. The concept of mass or density is discarded, the Equivalence Principle reflects only the equivalence of inertial effects and gravitation. An essential modification concerns the inertial terms which are written, like all the other terms of the equation of motion, in the form of a Helmholtz-Hodge decomposition [START_REF] Caltagirone | On Helmholtz-Hodge decomposition of inertia on a discrete local frame of reference[END_REF]. The equation of motion conceived both in discrete and continuous form is written:

∂v ∂t = -∇ Å φ o + 1 2 |v| 2 -c 2 l dt ∇ • v ã + ∇ ⊗ Å ψ o + 1 2 |v| 2 n -c 2 t dt ∇ × v ã + h s (3)
where h s = ∇φ s -∇ ⊗ ψ s is a source term also expressed as a Helmholtz-Hodge decomposition. Solving the equation ( 3) allows us to find the current solution at time t o + dt, v from its value v o and the equilibrium potentials φ o and ψ o at time t o . The divergence and curl of v allow the updating of the potentials in an explicit form:

updates            α l φ o -c 2 l dt ∇ • (v -v o ) -→ φ o α t ψ o -c 2 t dt ∇ × (v -v o ) -→ ψ o u + v dt -→ u o (4)
The symbol -→ corresponds to this update of ∇

• (v -v o ) and of ∇ × (v -v o ).
It is finally possible to update the displacement u. The parameters α l and α t are attenuation factors for longitudinal and transverse waves respectively. For example, if the longitudinal waves are not attenuated, α l = 1 and, similarly, if the transverse waves are completely dissipated, α t = 0. Mechanical dissipation in the form of heat by attenuation of compressional and shear waves completes the conservation of total energy fixed by the equation of motion. The dissipated energy between times t o and t o + dt is therefore

(1 -α l ) φ o + (1 -α t ) ψ o .
In the case where part of the compressional or shear energy is lost by viscous dissipation, the equation ( 3) must retain the full energy between the times t o and t o +dt where the velocity evolves between the equilibrium state v o and the current state v hence the notation (vv o ) in updates. For example, a rigid rotating flow, where the velocity of this incompressible flow is v r = Ω × r, does not correspond to energy dissipation since its primary curl is equal to ∇ × v r = cte; this curl cannot be accumulated in time within the vector potential.

The law of motion ( 3) is a wave equation that applies to the motion of fluids or solids as well as to waves of any kind, such as swell, acoustic waves or light. It is necessary to distinguish and separate once and for all the motions of the media to which they apply:

• the motion can be (i) non-zero rotational and of non-zero divergence, (ii) irrotational and compressible, (iii) non-zero rotational and solenoidal, and (iv) curl-free and divergence-free; there is, a priori, no other choice;

• the propagation medium can be a fluid, a solid or a vacuum; the equation (3) contains no constitutive law, state law or rheological law. The medium is defined only by its propagation properties, the longitudinal c l and transverse c t celerities and the respective attenuations of these waves, α l and α t .

This second part of the classification of a medium is not of particular interest, because fluids can be any, Newtonian, viscoelastic, viscoplastic, with threshold, time-dependent properties, etc. Equations that incorporate the constitutive laws of media or the characteristics of these media are confusing and lead to a misinterpretation of physical phenomena. This is the case of some equations which introduces physical characteristics such as density and dynamic viscosity which mix the two notions presented by considering, for example, that a constant density leads to an incompressible flow, or that a zero viscosity leads to an irrotational motion. The hyperbolic law of motion based on waves must contain only longitudinal and transverse wave propagation properties.

The form of the law (3) suggests that the motion can be decomposed into a divergence-free part and a curl-free part by separating the two contributions to the acceleration, γ = γ φ + γ ψ where γ φ = -∇φ and γ ψ = ∇ ⊗ ψ. The velocity then naturally decomposes into an equivalent form, v = v φ + v ψ but the scalar and vector potentials of the velocity will not be useful or sought. This splitting of the equation of motion is however not completely possible because of the nonlinear terms; indeed, the inertial potential |v| 2 /2 can be written as:

1 2 |v| 2 = 1 2 |v φ | 2 + 1 2 |v ψ | 2 + v φ • v ψ (5)
The quadratic terms |v φ | 2 /2 and |v ψ | 2 /2 remain in the strategy of decoupling compression and rotation effects and can be associated to Bernoulli potentials

φ o B = φ o + |v ψ | 2 /2 and ψ o B = ψ + |v ψ | 2 /2.
The equation (3) then reads as follows:

∂v ∂t + ∇ (v φ • v ψ ) -∇ ⊗ (v φ • v ψ n) = -∇ Å φ o + 1 2 |v φ | 2 -c 2 l dt ∇ • v φ ã +∇ ⊗ Å ψ o + 1 2 |v ψ | 2 n -c 2 t dt ∇ × v ψ ã + ∇φ s -∇ ⊗ ψ s (6)
On the other hand, the terms (v φ • v ψ ) prevent any possibility of separating the irrotational and solenoidal components, unless one of the components is zero. In the general case, the compressional and rotational effects are intertwined and the equation ( 6) must be solved without simplification. An attempt by D.D. Joseph [START_REF] Joseph | Helmholtz decomposition coupling rotational to irrotational flow of a viscous fluid[END_REF] to decouple the irrotational terms from the solenoidal terms of the Navier-Stokes equation shows that this Helmholtz-Hodge decomposition is legitimate to obtain a solution suitable for the type of solution sought; however, the Navier-Stokes equation has characteristics that prevent this operation.

Let us now consider the cases where the equation ( 6) can be disassembled; there are only two possibilities, (i) the one corresponding to a velocity curl-free, ∇ × v = 0 where the velocity field reduces to its only irrotational component v = v φ (ii) the one where it is the divergence of the velocity which is null ∇ • v = 0 and where the velocity field is solenoidal v = v ψ . The case where the solution is harmonic, i.e. without curl or divergence, does not exist. As [START_REF] Caltagirone | Extension of galilean invariance to uniform motions for a relativistic equation of fluid flows[END_REF] shows, these solutions with zero or constant divergence or zero or constant curl are structurally discarded from the solutions of the equation [START_REF] Peyret | Computational Methods for Fluid Flow[END_REF].

When ∇ × v = 0, the equation of motion can be reduced to finding its unique irrotational solution, v φ , it is written:

       ∂v φ ∂t = -∇ Å φ o + 1 2 |v φ | 2 -dt c 2 l ∇ • v φ ã -∇φ s α l φ o -dt c 2 l ∇ • v φ -v o φ -→ φ o (7)
This equation ( 7) is equivalent to the Euler equation of motion where the fluid is ideal and the flow compressible where its solution is of the form (v, φ). However, the differences are important: (i) the presented formulation has no constitutive law, the state law is contained in the scalar potential φ = p/ρ, (ii) the thermodynamic evolution is integrated in the longitudinal celerity c l which can reflect an isothermal or isentropic evolution, (iii) it does not require an associated energy equation because the equation ( 7) already expresses the conservation of the compression energy and where φ represents the mechanical energy of compression per unit of mass. Of course, the fluid considered cannot be viscous nor conduct heat, it is necessarily an ideal fluid. This formalism can potentially represent compressible flows without shocks or in presence of shocks. In the latter case, a term representing the discontinuity of the form ∇(ξ φ), where ξ is a phase function, must be added to the equation of motion. Figure 1(a) schematizes the associated primary structure where the flow is basically directional although the overall domain is threedimensional. The velocity divergence ∇ • v φ reflects only the cause-and-effect interactions from one Γ segment to the next, the interlocking of direct and induced currents has been broken and propagates only through the vertices. When ∇ • v = 0, the solenoidal solution of the equation of motion can be found directly from component v ψ :

       ∂v ψ ∂t = ∇ ⊗ Å ψ o + 1 2 |v ψ | 2 n -dt c 2 t ∇ × v ψ ã -∇ ⊗ ψ s α t ψ o -dt c 2 t ∇ × v ψ -v o ψ -→ ψ o (8)
where α t is the attenuation factor of transverse waves. The formulation can be used without a constitutive law, it is dt c 2 t = ν, the kinematic viscosity which defines the viscous interactions within the fluid; it can be any function of another variable, for example the temperature or directly the shear energy ψ o . It can depend on the direction considered since this transverse propagation term allows to take into account anisotropic and polarizable media. In the case of an isotropic and homogeneous medium such as a Newtonian viscous fluid, the attenuation of transverse waves is fully and α t = 0, the instantaneous shear stress is equal to

ψ o = ν ∇ × Ä v ψ -v o ψ ä
, the energy lost by viscous dissipation and transformed into heat is thus equal to

(1 -α t ) ψ o .
The solutions of the equation ( 3) are also those of the Navier-Stokes equation when ∇ • V = 0 is imposed with sufficient accuracy, many simulations show clearly [START_REF] Caltagirone | Application of discrete mechanics model to jump conditions in two-phase flows[END_REF]. In this case, the mechanical equilibrium is ensured by the term ∇φ o . What changes in the equation ( 8) is that the equilibrium is assured by the term ∇ ⊗ ψ o but the solution of the problem is strictly the same. Thus the relation -∇φ o = ∇ ⊗ ψ o makes it possible to adapt the equilibrium term to the problem posed, a flow with zero divergence without the presence of the scalar potential (pressure) reserved for compressible flows. The boundary conditions of the problem posed are not formulated in the same way as in continuum mechanics, in the form of partial derivatives on the edge, but in terms of accelerations directly integrated in the equation of motion [START_REF] Caltagirone | Generalized local boundary conditions in terms of accelerations on an edge-based formulation[END_REF]. They therefore be absolutely arbitrary since they correspond to fluid flows or to localized shear stresses for each part of the domain boundary.

The law of motion (3) has remarkable properties related to its structure which is a decomposition of the intrinsic acceleration. It is naturally relativistic [START_REF] Caltagirone | Extension of galilean invariance to uniform motions for a relativistic equation of fluid flows[END_REF] because it is a non-linear wave equation; it satisfies Lorentz invariance without using the factor γ = 1/ » 1v 2 /c 2 0 , c 0 being here the celerity of light. In classical applications, the velocity is negligible compared to c 0 , but the theory of relativity cannot be ignored in the derivation of any new physical law.

Another event that marked the science of the 20th century is Noether's theorem [START_REF] Kosmann-Schwarzbach | Noether Theorems. Invariance and Conservations Laws[END_REF] which states that if a physical law can be expressed in the form of a Lagrangian or Hamiltonian, then it has the properties of invariance and conservation of some quantities. The discrete law is formed by two Lagrangians associated respectively with longitudinal and transverse effects; each of these Lagrangians has two actions represented by the potential and kinetic energies. Thus, although the equation ( 3) is formulated in terms of accelerations, the time invariance is related to the conservation of the total energy, the translational invariance to the conservation of the compressive energy and the rotational invariance to the conservation of the angular acceleration momentum in continuum mechanics).

Primal and dual local form of Bernoulli's theorem

The two-term decomposed form, without divergence and without curl, of the equation of motion suggests a local form on a segment Γ equivalent to Bernoulli's theorem for an ideal fluid. In discrete mechanics, the local form of the theorem is directly deduced from the compression energy per unit mass defined by the gradient of the scalar potential φ of the acceleration. Indeed, the integral along the segment Γ of the primal structure of the acceleration represents the compression energy:

ˆΓ γ φ • t dl = -ˆΓ ∇φ • t dl (9)
When the length of the segment tends to zero, we obtain the equivalent local form of Bernoulli's theorem:

∂v φ ∂t = -∇ Å φ o + 1 2 |v φ | 2 -dt c 2 l ∇ • v φ -φ s ã along the path Γ and in Ω ( 10
)
where φ s is a scalar potential, for example φ s = g z for gravity. When the motion is irrotational, it is possible to traverse the whole Ω domain from one vertex of the primary structure to the other by a discrete integration on the Γ segment. The path once defined allows to cover the whole physical domain. If the first two terms are identical, the differences between this form [START_REF] Caltagirone | Generalized local boundary conditions in terms of accelerations on an edge-based formulation[END_REF] and the classical form of Bernoulli's theorem are multiple [START_REF] Landau | Fluid Mechanics[END_REF][START_REF] Truesdell | A First Course in Rational Continuum Mechanics[END_REF][START_REF] Truesdell | Daniel Bernoulli's Hydrodynamica[END_REF], (i) the flow can be unsteady and (ii) the flow can be compressible with a term -dt c 2 l ∇ • v φ which is expressed as a function of the velocity itself. The statement of Bernoulli's theorem for irrotational flows applies equally well to a Lagrangian and an Eulerian description:

Primal Bernoulli theorem: If, at a given time, an irrotational flow of an ideal fluid is compressible, inertial and derived from a scalar potential, it will remain so during its evolution in time.

An equivalent statement can be deduced to represent the evolution of the compression energy unit mass on the segment Γ:

φ o + 1 2 |v φ | 2 -dt c 2 l ∇ • v φ -φ s = cte along the path Γ ( 11
)
where the term dt c 2 l ∇ • v φ , also expressed on the segment, corresponds to the transformation of kinetic energy into potential energy symbolized by φ o .

Bernoulli's law ( 10) is naturally relativistic [START_REF] Caltagirone | Extension of galilean invariance to uniform motions for a relativistic equation of fluid flows[END_REF] without applying a priori the Lorentz transformation. The most important consequence concerns the limit velocity that a fluid can reach by following a rectilinear trajectory where it cannot exceed the celerity of the medium, v < c 0 for light in vacuum or v < c l for fluids. The most significant example is that of the Sod's tube [START_REF] Sedov | Similarity and Dimensional Methods in Mechanics[END_REF][START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF] where the velocity of the gas contained in a channel cannot exceed the local celerity of sound; a shock wave materializes the propagation of the wave front in the pipe.

In the presented formalism, φ and ψ are the acceleration potentials, two energies per unit mass. In the discrete incompressible formulation, it is the energy ψ associated with the induced effects that must be considered and not the velocity:

ˆΓ γ ψ • t dl = ˆΓ ∇ ⊗ ψ • t dl ( 12
)
or a local form:

∂v ψ ∂t = ∇ ⊗ Å ψ o + 1 2 |v ψ | 2 n -dt c 2 t ∇ × v ψ -ψ s ã along the path Γ and in Ω (13)
This relation ( 13) applies with a minimum of assumptions; the flow can be unsteady, inertial, no-dissipative and derive from a vector source potential ψ s . An equivalent statement can be deduced to represent the evolution of the rotation energy per unit mass on closed contour ∆:

ψ o + 1 2 |v ψ | 2 n -dt c 2 t ∇ × v ψ -ψ s = cte along the closed path ∆ ( 14
)
where the term dtc 2 t ∇×v ψ , also expressed on the contour ∆, corresponds to the transformation of transverse kinetic energy in potential energy or its inverse. The statement of Bernoulli's theorem in its dual form is therefore the following:

Dual Bernoulli theorem: If, at a given time, a solenoidal flow of an incompressible fluid is inertial and derived from a vector potential, it will remain so during its evolution in time.

The equations ( 10) and ( 13) have a perfectly symmetric structure. The conditions of application of the equation ( 10) are simple, it is sufficient that ∇ × v = 0 for the motion to be irrotational. For equation [START_REF] Truesdell | A First Course in Rational Continuum Mechanics[END_REF], it is sufficient that ∇ • v = 0 for the flow to be solenoidal. This result is due to the absence of harmonic motion, both curl-free and divergence-free [START_REF] Caltagirone | Extension of galilean invariance to uniform motions for a relativistic equation of fluid flows[END_REF].

By returning to the complete form of the law of motion, it is possible to represent the entanglement of the compressional and rotational effects on a Γ trajectory; this form extends, as for the two Bernoulli theorems, to the whole physical domain. Figure 2 schematizes the entanglement of these irrotational and solenoidal terms. This diagram illustrates very well the J.C. Maxwell's inspiration [START_REF] Maxwell | A dynamical theory of the electromagnetic field[END_REF] to couple the electrical and magnetic effects in the same system of equations. The link between Bernoulli's primal and dual theorems can be interpreted in the light of Maxwell's remarkable idea of representing direct and induced currents on the same support Γ.

An equation of motion for incompressible flows of viscous fluid

The Navier-Stokes equation is always accompanied by a conservation of mass constraint, the continuity equation in the general case or ∇ • V = 0 the incompressible case. If the divergence were zero at any time, it would be impossible to obtain the pressure field; it is therefore the residual divergence ∇ • V ≈ ε that allows us to find the value of the pressure in order to ensure mechanical equilibrium. Similarly, in discrete mechanics the scalar potential of the acceleration φ o = p/ρ is obtained by the accumulation process of the equation (3) from ∇ • v. In both cases, the equilibrium of the equation of motion associated with a zero or smallest possible residue requires the knowledge of ∇ • v.

The proposed formulation consists in noticing that the scalar potential φ (or the pressure p = ρ φ) can be replaced by a vector potential of the acceleration ψ. Any acceleration can be written as a gradient of a scalar function or as a dual curl of a vector function:

-∇φ = ∇ ⊗ ψ (15)
where φ and ψ are two functions generally of different variables. The equation of motion in its incompressible acceptation is then:

       ∂v ψ ∂t = ∇ ⊗ Å ψ o + 1 2 |v ψ | 2 n -ν ∇ × v ψ ã ψ o -ν ∇ × v ψ -v o ψ -v o ψ -→ ψ o (16)
The uniform rotation flow must be filtered beforehand into the energy accumulation so that

ν ∇ × Ä v ψ -v o ψ ä
represents only shear and not full rotation. This distinction for Newtonian fluids is required by the complete dissipation of the shear energy between t o and t o + dt. The partial derivative in time can be evaluated at any order, but assume here that it is represented by a first-order Euler-type scheme; the current velocity v ψ is then written as a function of that obtained at the previous instant v o ψ :

v ψ = v o ψ dt ∇ ⊗ Å ψ o + 1 2 |v ψ | 2 n -ν ∇ × v ψ ã (17)
Given the identity ∇ • (∇ ⊗ ψ) = 0, the velocity v ψ remains at zero divergence. However, it is necessary that ∇ • v o ψ = 0 is initially at zero divergence, including the boundary conditions attached to it. The simulations performed with the equation ( 16) show indeed that the divergence remains at the level of the accuracy of the machine throughout the computation. This equation can be treated explicitly, in its form [START_REF] Maxwell | A dynamical theory of the electromagnetic field[END_REF] or implicitly by semi-implicating the nonlinear term 3 represents schematically the process of solving in time the equation [START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF]. From the value of the vector potential ψ o , for example zero at the initial instant, the implicit resolution of the equation provides the velocity v ψ ; from the curl of this one the potential ψ o is explicitly re-actualized and a new cycle allows to progress in time. The interest of the presented formulation is obvious, it allows to suppress the projection step on a field of zero divergence of the formulations where the pressure is used as for the Navier-Stokes equation. The resolution of this equation of motion requires the application of a specific method to reduce or even cancel the divergence of the prediction step. This is the case of the family of methods based on a Poisson equation on the pressure [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. Implicit or semi-implicit methods such as the augmented Lagrange method are also efficient but more expensive in terms of computation time [START_REF] Fortin | Méthodes de lagrangien augmenté, Application à la résolution numérique de problèmes aux limites[END_REF]. In any case, the method consisting in imposing a zero or limited divergence represents a very important part of the total time of solution of the Navier-Stokes equation.

|v ψ | 2 = v ψ • v o ψ . Figure
Equation ( 17) is remarkably fast to solve because the field v ψ must be searched in an already divergence-free space where the property ∇ • (∇ ⊗ ψ) = 0 is satisfied structurally. The equivalent primary operator ∇(∇ • v) = -∇(∇ • v) has been used to project a velocity field v ψ = v + v onto a zero divergence field where v satisfies well ∇ • v = -∇ • v. This operator developed by several authors from 1999 onwards [START_REF] Caltagirone | A vector projection method for solving the Navier-Stokes equations[END_REF] for the method of vector projection shows that, if the second member is of the same nature as the first, then the solution is obtained by solving the linear system in a few iterations of a conjugate gradient method, independently of the number of unknowns. This exceptional convergence is explained by the algebraic properties of the linear system and of the spectrum of the discrete operator.

The validation of this approach will require a long validation work in different situations. In order to show the interest of the formulation, only two simple cases are analyzed, the plane Poiseuille flow and a flow in a vertical channel subjected to a temperature difference.

Two elementary examples of solenoidal flows

The two proposed examples cover only a small part of the potential of the presented formulation. The corresponding computer simulations have been well performed using the equation ( 16) and show the expected results, (i) the solution of the Poiseuille problem is obtained exactly (to machine precision) since the method is second order in space (ii) the solution of the heated channel is second order since the solution corresponds to a polynomial of degree three. Given the theoretical nature of the background, only explanations of the physical model are developed.

Planar Poiseuille flow

The Poiseuille stationary flow of a viscous fluid in a planar channel is emblematic of those in fluid mechanics for incompressible flows. Let us consider a channel of height h = ±1 in which a fluid of kinematic viscosity ν flows under the effect of a constraint, a pressure gradient or an imposed velocity, here equal to V 0 = 1, the mean velocity in the channel.

In continuum mechanics, the solution is obtained from the Navier-Stokes equation where the inertia term takes several forms, including V

• ∇V = ∇(V 2 /2) -V × ∇ × V: ρ ∇ Å 1 2 V 2 ã -ρ V × ∇ × V = -∇p + ∇ × (µ ∇ × V) (18)
solution is represented by the velocity, V • e x = u(y) = 3/2 1 -y 2 and by pressure p(x) = -3 µ x to a constant or by p B (x, y) = p(x) + 1/2 ρ V 2 = -3 µ x + 9/8 ρ (1 -y 2 ) 2 , Bernoulli pressure. The projection onto the (x, y) components reads as follows:

9 2 ρ (1 -y 2 )y e y - 9 2 ρ (1 -y 2 )y e y = - ∂p(x) ∂x -3 µ e x (19)
Thus, the pressure gradient balances the viscous terms and the inertia terms cancel, the Poiseuille flow is non-inertial but only by compensation. Note that in this case the Lamb vector of the form -V×∇×V is a gradient and not a curl. In the general case, its form is indeterminate.

In the discrete formalism the equation of motion reads:

∇ Å |v| 2 2 ã -∇ ⊗ Å |v| 2 2 n ã = ∇ ⊗ ψ o -∇ ⊗ (ν ∇ × v) (20)
The solution corresponds to the equilibrium of the Poiseuille flow:

9 2 (1 -y 2 )y t - 9 2 (1 -y 2 )y t = -∇ ⊗ ψ o (y) -3 ν t (21)
or a vector potential equal to ψ o (y) = 3 ν y n or ψ o B (y) = 3 ν y n + 9/8 (1 -y 2 ) 2 n the Bernoulli vector potential.

These solutions of the continuum mechanics and the discrete formulation lead of course to the same result with scalar and vector potentials of different functions:

-∇φ o (x) = ∇ ⊗ ψ o (y) = 3 ν e x (22)
The functions -∇φ o (x) and ∇⊗ψ o (y) of two different variables are both equal to a constant. This is a special case because, in general, these two functions are arbitrary but the equality -∇φ o = ∇ ⊗ ψ o can always be used to select the appropriate potential.

Vertical channel differentially heated

Consider a vertical planar channel whose no-slip boundaries define a fluid cavity where x ∈ [-h, h] and the direction y is undefined; the viscous, expandable fluid is heated differentially with a temperature difference ∆T between the vertical planes. The gravitational force acting on each portion of the fluid is due to the dilatability of the fluid with temperature. Using the Boussinesq approximation, this force per unit mass is written:

ρ ρ 0 g = - Å g δ 2 x h ã e y = -∇φ h + ∇ ⊗ ψ h (23)
where φ h and ψ h are the Archimedean potentials and δ = β∆T where β is the thermal expansion coefficient.

The Navier-Stokes equation or the presented model lead of course to the same solution,

v(x) = g δ h 2 12 ν Å x h - x 3 h 3 ã e y , (24)
where, for a low enough Rayleigh number, the inertia has been neglected. The scalar and vector potentials, defined with respect to a constant, can be written in terms of the coordinates (x, y) as:

         φ h (x, y) = g δ 2
x y h

ψ h (x) = - g δ 4 x 2 h e z (25)
In the case where the scalar potential is chosen to balance the equation of motion, it is expressed as a function of both spatial coordinates. When the formulation is based on the vector potential, it depends on only one coordinate. This difference reflects the consistency of the formulation (v, ψ), as the physical problem posed does not depend on the vertical coordinate at all. It is to be noted that the equation ( 16) provides, at the same time, v and ψ with ψ = ν ∇×v.

These two simple examples do not validate the physical model presented, but they allow us to understand why it is more appropriate to use a velocity-vector potential formulation of the acceleration (v, ψ) an incompressible flow.

Discussion

The consequences of separating acceleration and velocity into solenoidal and irrotational components cannot be drawn immediately. For example, the impact on the boundary conditions would deserve further investigation. However, the superposition of direct and induced currents on the same medium suggested by Maxwell's theory leads, for viscous fluids, to a dual Bernoulli theorem where the sum of the divergence-free terms is equal to a constant along a closed path. The equation of motion specific to solenoidal flows benefits from the form of this theorem where its solution remains divergence-free if it is divergence-free at a given time. From a practical point of view, the solution of this equation does not require to project the current solution on a field of zero divergence if the solution at the instant preceding the instant t o is already of zero divergence. The difficulty of transforming a gradient field into an equivalent curl field is real in fluid mechanics equations and in particular for the Navier-Stokes equation. The proposed formulation suppresses the projection step in the case where only the solenoidal solution is sought.

A real flow can be both to non-zero divergence and non-zero curl, but it is not necessarily turbulent; it can also be inertial and non-turbulent. One of the questions that can legitimately be asked is: can a strictly incompressible flow be turbulent? The answer is not found among the many simulations of turbulence performed from the Navier-Stokes equation because, when solving it, at certain stages, the flow is not at divergence-free. In general, if the pressure can be calculated, it is because the divergence was not zero at a given time. The constraint ∇ • V = 0 cannot be imposed a priori in the Navier-Stokes equation. The presented formulation is strictly solenoidal and the problem remains open.

Under what conditions a compressible, viscous, inertial flow become turbulent? The turbulence statistics of computationally simulated flows will not tell us either. The respective roles of rotational and compressive effects are not yet known with precision. In particular, the effects of acoustic propagation of longitudinal waves are largely underestimated at present. The separation of divergence-free and curl-free terms in the velocity field is a new way to study physical turbulence. Solving the equation of motion in two components v φ and v ψ will probably allow a finer analysis of the turbulent character of a flow. The inertial terms in v 2 φ and v 2 ψ correspond respectively to the pressure and the shear energy which do not interact with each other and can simply be integrated in the potentials to form the Bernoulli potentials. On the other hand, the inertial terms (v φ • v ψ ) alone allow the mutual transfer of compressional and rotational energies over time. It seems possible that the product (v φ •v ψ ) represents the cornerstone of understanding the physics of turbulent flows.

Declaration of competing interest

There are no conflict of interest in this work.

Figure 1 .

 1 Figure 1. Primal (a) and dual (b) structures defined respectively by the contours Γ and ∆. The primal contour Γ consists of segments of length dh = [a, b] oriented by the unit vector t; the scalar potential φ is positioned on the vertices of the segments. The closed dual contour ∆ is oriented by the unit vector n where the condition n • t = 0 is realized by construction; the vector potential ψ is carried by the vector n. The primal contour Γ can be open (a) or closed (b). The acceleration and velocities are expressed on the segment Γ.

Figure 2 .

 2 Figure 2. Primal (a) and dual (b) structures defined respectively by contours Γ and ∆ mimicking the dynamic entanglement between direct and induced currents; each segment Γ of length dh = [a, b] is oriented by a unit vector t. The normals to the primal facets are also oriented along n with n • t = 0. The scalar potential φ is defined on each vertex of this primitive structure and the vector potential ψ is carried by unit vector n. The acceleration and velocities are expressed on the segment Γ.

Figure 3 .

 3 Figure 3. Diagram of the equation solving process[START_REF] Sod | A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws[END_REF]


Luc Absil

Yann Balland

Franck Pereira
email: franck.pereira@obspm.fr

Dos Santos

Long-range temperature-controlled transport of ultra-cold atoms with an accelerated lattice

We report our method for transporting ultracold atoms over macroscopic distances and trapping them back in a vertical mixed trap, consisting of the superposition of a vertical lattice and a transverse confinement beam. The transport is performed with Bloch oscillations allowing us to move up to 25% of a sub-micro-Kelvin atomic cloud on a distance of the order of 30cm, without excessive heating and with a good control of its final position. The efficiency is lowered to about 10% after trapping them back in the vertical mixed trap during extended times.

I. INTRODUCTION

The transport of clouds of cold atoms or Bose-Einstein condensates (BEC) from their production area to other areas of experimental interest has drawn much attention since first implementations in the 90s [1]. Various approaches have been considered, using magnetic or optical fields, depending on the desired specifications. Magnetic traps have demonstrated to be a suitable solution for long-range displacement of tens of centimeters, but at the cost of heating up the sample [2]. On the other hand, trapping the atomic cloud in optical lattices has proven to be a solution of choice for precise displacement of single atoms over distances of the order of the centimeter [3] [4], for adiabatic centimeter-long transport of ultracold atoms [5]. Other solutions involving the use of Bessel beams [6], of a focus-tunable moir lens [7], other optical systems of tunable-zoom [8] or optical tweezers [9] have also shown their interest to transport atomic cloud either on short distances with reasonably small heating or on long distances with more consequent warm-up of the atomic cloud.

In the experiment reported here, we aim at performing ultra-sensitive short range forces measurements, with a quantum sensor based on 87 Rb atoms trapped in an optical trap, consisting of a vertical shallow lattice λ V erdi = 532 nm overlapped with a transverse confinement potential λ IR = 1064 nm. In such a system, energy differences between adjacent lattice sites can be probed with a Ramsey-Raman interferometer [10] providing us a local measurement of the vertical force exerted on the atoms. The goal is to probe micrometer-range forces occurring between our atomic cloud and the surface of a mirror. Previous attempts at similar measures have shown the importance of preparing the sample away from the surface of interest in order to minimize any undesired adsorption of the atoms onto the surface [11] [12]. Furthermore, the spatial resolution of the measurements is limited by the final size of the cloud, in the first place in the vertical direction, perpendicular to the surface of the mirror, of the order of a few micrometers [10]. Secondly, the transverse elongation of the cloud leads to an averaging of the roughness of the surface.

We present here the outcomes of our experimental setup which enables the transportation of up to 25% of a sub-micro-Kelvin atomic cloud on a distance of the order of 30 cm, without inducing excessive heating and while maintaining control over its size and final position.

We first provide a general overview of the experimental setup and of the main sources of loss expected, before reporting our results on the transport efficiency and the control of the atoms position.

II. EXPERIMENTAL SETUP

Figure 1-a summarizes the main steps of the atom preparation. A two-dimensional magneto-optical-trap (2D-MOT) loads atoms in a three-dimensional magneto-optical trap (3D-MOT), which in turn loads a crossed dipole trap where evaporative cooling is performed down to a temperature of 300 nK and an atomic density of the order of 10 12 atoms/cm 3 .

Then, atoms are launched upward and transported over 30 cm, up to the vicinity of the surface of the lattice mirror, where measurements are to be performed.

The elevator consists of the overlapping of two counterpropagating beams "Bloch 1 " and "Bloch 2 " whose frequency differences and intensities are tuned over time. Atoms are thus loaded in an accelerated or decelerated lattice of adjustable depth.

A. Setup

The Bloch beams are derived from a single external-cavity diode laser, which is split in two. Each beam is sent to a double pass acousto-optic modulator (AOM) before injecting a 2W tapered amplifier. Each amplified beam then passes through an other single pass AOM and then injects a polarization maintaining fiber. One fiber directs the first Bloch beam to the bottom of the experimental setup, where it enters the vacuum chamber, while the other fiber bring the second Bloch beam to the top. They are finally overlapped so as to counterpropagate along the vertical direction. This setup allows for modulating the frequency difference between the beams during the transport with the AOMs before the TAs, and their powers with the AOMs placed after the TAs.

B. Atoms in moving lattices

The periodic potential U (z) resulting from the interference of the two counter-propagating beams "Bloch 1,2 " of intensity I 1,2 is proportional to the local intensity experienced by the atoms I(z) = 2I 0 (1 + m cos 2kz), where I 0 is the average intensity and m the ratio of the geometric and arithmetic means of the intensities of each beam :

U (z) = U 0 2 (1 + m cos 2kz) (1)
For a detuning ∆ of the order of hundreds of GHz with respect to the D2 line of 87 Rb, the depth is U 0 ∼ hΓ 2 I 0 /(3∆I s ), where Γ and I s are the natural linewidth and the saturation intensity of the transition, and k is the modulus of the wavevector of the laser. By imposing a frequency difference ∆ω between the two Bloch beams, we then move the lattice potential with a velocity v latt = ∆ω/2k. This allows us to coherently accelerate the atoms through Bloch oscillations [13], each oscillation imparting a 2hk momentum transfer to the atoms.

For the measurements presented here, the detuning ∆ was set to -250 GHz. For maximum powers of order of 50 mW in the two beams, with identical waists of 400 µm, we calculate a maximum depth of the order of 150E r . This is in reasonable agreement with a determination of the depth based on Raman-Nath diffraction, where from the measurement of the populations in the diffracted orders as a function of duration of a square pulse performed before the transport, we deduced a maximum depth of about 110E r .

Figure 2 illustrates the typical sequence of transport. The atoms are first adiabatically loaded into the lattice, by increasing the depth from 0 to U acc in 600 µs while ramping the frequency difference between the two beams by 25 MHz/s. This corresponds to accelerating the lattice by g.

Then, the lattice is accelerated upwards over 6 ms at a fixed lattice depth U acc . The chirp of the frequency difference between the beams corresponds to N acc Bloch oscillations. After a certain transport time at constant velocity, where the lattice potential may be reduced, the atoms are decelerated at a depth U dec with N dec Bloch oscillations. The atoms are then adiabatically released by reducing the lattice potential down to zero in 600 µs. They are then trapped in the mixed trap consisting of the superposition of the green lattice and the IR transverse confinement beam illustrated in figure 1, where force measurements are to be performed.

Between the acceleration and deceleration, we keep the lattice at a depth U f f in order to keep the atoms vertically and transversally confined. Additional transverse confinement is also provided by the IR laser at 1064 nm whose waist, of the order of 200 µm is placed at the target position of the atoms, near the surface of the lattice mirror.

C. Sources of loss

Common sources of loss during the transport with Bloch oscillations are spontaneous emission and non-adiabatic losses. Spontaneous emission occurs at a rate of order of P sp ∼ U 0 Γ/h∆ and leads to losses that increase with the duration of the transport. Note that in our case, the tapered amplifiers emit in addition amplified spontaneous emission near resonance, which we efficiently filter using narrow bandwidth (of order of ∼ 150 GHz) interference laser line filters. This in turn forces us to use large enough detunings.

For their efficient transport with Bloch oscillations, atoms have to be initially loaded and then to remain in the lower band of the moving lattice. The adiabaticity of the acceleration process restricts accelerations to below a critical value, which, in the weak binding limit, can be deduced from the Landau-Zener formula [13]. The larger the depth, the larger this critical acceleration, and generally losses decrease with decreasing accelerations, and thus longer transport durations, by contrast with losses by spontaneous emission. This calls for a compromize between these two sources of losses. Note that for depths below 20

E r , resonances between the energy bands also reduces the efficiency in certain acceleration ranges.

Finally, specific to our optical setup, a last source of loss is due to small fractions of the Bloch 1 beam that retroreflect on the green lattice and Raman mirrors. These create additional static lattices, in which a fraction of the atoms, can be trapped as displayed in subfigure 1-α. To minimize them, on one hand, the dichroic character of the green lattice mirror yields a low reflectivity of about 0.6% at 780 nm, and on the other hand, the polarization beam splitter placed on the top should ideally prevents the Bloch 1 beam from reflecting onto the Raman mirror. In practice, a residual reflection is minimized by carefully adjusting the polarization of the beam.

D. Detection

Two different detection schemes are used to evaluate the performance of our transport system. We use in-situ absorption imaging with a PCO PixelFly CCD camera in the bottom chamber, and both absorption imaging or fluorescence detection with a deeply cooled EM-CCD camera in the top chamber. This allows us to measure the number of atoms in the |F = 2 state and their temperature in both chambers. In-situ fluorescence imaging in the top chamber allows to measure with a better atom detectivity but without much spatial resolution the number of atoms in each of the |F = 1 and |F = 2 states.

III. RESULTS

To optimize the launch efficiency of the atoms, we started by imaging in the bottom chamber their position 3 ms after their acceleration. Subfigure 1-α left displays an absorption imaging picture featuring three different groups of atoms: (3) launched, (2) trapped in a parasitic lattice due reflection of the "Bloch 1" beam on the Raman mirror, and

(1) untrapped and in free fall. By adjusting the polarization of the "Bloch 1" beam, we eliminated the loss of atoms in both groups (1) and ( 2), which results in a single group of atoms as displayed at the right. We then increased the number of Bloch oscillations until the atoms could reach the top chamber, where we obtained a first signal (similar to the one in subfigure 1-β). We then optimized the durations of the different transport phases, keeping the transport distance constant. This initial optimization resulted in durations of 6 ms for both acceleration and deceleration phases, with 250 Bloch oscillations at the acceleration as displayed in figure 2. We then optimized the other parameters, such as the power of the "Bloch" beams or the temperature of the atoms before the transport. We finally varied back the duration and number of Bloch oscillations without much improvement.

We illustrate Fig. 3 the influence of the initial temperature of the atoms both on the number and the temperature of the atoms after their transport over the 30cm distance and after their trapping in the mixed trap described above. We find that the efficiency of both processes reduces when increasing the initial temperature. We attribute this for one to the velocity selectivity of the acceleration process, that works best for atoms moving at subrecoil velocities in the frame of the lattice, and for the other to the shallow depth of the mixed trap. Note that the efficiency is not the only relevant parameter when it comes to selecting the most appropriate parameters for our transport system. It can be of interest for instance to work either with high initial temperatures, which, despite the lower efficiency, leads indeed to larger number of transported atoms, or with lower initial temperature, which leads to smaller sizes at the end. These relevant parameters actually depend on the details of the evaporative cooling stage, where temperature, number of atoms and initial size are correlated.

The influence of the powers in the Bloch beams during the transport is illustrated in figure 4-a. The larger the powers, the higher the depth, and the more efficient is the transport.

However, for an efficient subsequent loading in the mixed trap, the power in the "Bloch 1 " beam needs to be reduced, as shown in Figure 4-c, which displays the number of atoms in the mixed trap after a trapping time of 1s. We attribute this feature to the presence of residual parasitic static lattices, related to the reflection of the "Bloch 1 " beam on the Raman and green lattice mirrors.

Other parameters have also been varied without much improvement. We can for instance mention the influence of the depth used during the transport (U f f ). While decreasing the depth at this stage reduces in principle losses due to spontaneous emission, it comes at the price of a large decrease in the number of transported atoms, which we attribute to the reduction in the transverse confinement. Therefore this depth was kept constant, as high as during the acceleration and deceleration phases, for all the measurements presented here.

The position stability of the cloud after the transport is a key objective for our experimental setup. To check its evolution over time, we alternated absorption imaging measurements of the cloud before and after the transport. The corresponding Allan standard deviations of those two signals and of their difference are plotted in figure 5.

A larger measurement noise is observed for the position before the launch in the bottom chamber, which we attribute to the larger size of the cloud, whose center position has been measured after a 6 ms thermal expansion in free fall, and to the larger readout noise of the bottom camera compared to the top one.

On both positions of the cloud, we observe a bump on the Allan deviation, occurring around 200s. This perturbation is correlated with the temperature cycles of our air conditioning system and is actually eliminated when calculating the difference of the two signals,

since the Allan deviation of this difference averages as white noise down to about 100 nm at 1000s. This confirms that the travel distance is not affected by these thermal fluctuations, which are dominated by the fluctuations of the starting position of the atomic cloud, before the transport.

This correlation between the top and bottom positions decreases over longer time scales, of hours or days, as can already be observed in figure 5 for timescales of a few thousand seconds. This could actually be due to deformations in the mechanical structure supporting our two imaging systems, rather than variations in the actual travel distance.

Finally, the final vertical size of the cloud will determine the spatial resolution of our force sensor. But, given the current resolution of our imaging system of about 9µm, we are actually not able to precisely determine this size, nor its fluctuations. Nevertheless, since the atoms remain trapped at all times, we expect the size of the atomic cloud to be identical to the size at the end of the evaporation (ie of order of a few micrometers in rms width).

IV. CONCLUSION

We have investigated the efficiency of using Bloch oscillations to transport a cloud of ultracold atoms over a long-range distance of 30 cm, with an accurate control of their final position, size and temperature, and with the strong constraint of keeping them trapped at all times. We obtain a maximum efficiency of 25% for their transport, and in the end, we recapture about 10% of the initial number of atoms in the final trap, which combines a blue-detuned static lattice and a red-detuned progressive wave for transverse confinement.

Spontaneous emission induces in our setup significant losses and heating, which could be reduced by operating at larger detunings, and thus intensities. At the end of the sequence, we are left with a cold and well localized sample of atoms, which we will use to probe short range interaction forces between the atoms and the dielectric mirror that retroreflects the lattice beam.

1 FIG. 1 :FIG. 2 :

 112 FIG.1: Experimental setup of our experiment.Subfigures [α] and [β] provide examples of absorption images of the atomic cloud, respectively 3msafter the starting of the elevator and 30cm higher, at the end of the elevator, near the surface of interest. In [α] two situations are illustrated, either in the case of an imperfect (left) polarization, or in the case of an adjusted (right) polarization of the "Bloch 1 " beam, which emphasizes the losses due to the parasitic lattice due to the reflection of the "Bloch 1 " beam on the Raman mirror ( 2 ) and the losses due to the inadiabaticity of the launch ( 1 ).

FIG. 3 :

 3 FIG.3: Influence of the initial temperature of the cloud (at the end of the cooling processes and right before the launch) on (top) the number of atoms at the end of the transport and after 1s in the mixed trap and (bottom) the temperature in the vertical direction.

FIG. 4 :

 4 FIG. 4: Influence of the powers in the Bloch beams on the number of atoms and vertical temperature T y , either at the end of transport (a and b), or after their recapture in the mixed trap (during 1s) (c and d).

FIG. 5 :

 5 FIG. 5: Overlapping Allan standard deviation of the center position of the atomic cloud before and after the transport


Yanandlall Gopee
email: ygopee@laas.fr

Guy Camilleri
email: guy.camilleri@irit.fr

Margot Gaetani-Liseo

Xavier Roboam
email: xavier.roboam@laplace.univ-tlse.fr

Anne Blavette
email: anne.blavette@ens-rennes.fr

Corinne Alonso
email: alonsoc@laas.fr

Energy Management System for a Low Voltage Direct Current Microgrid: Modeling and experimental validation

Keywords: Microgrid, energy management system, optimization, linear programming, battery degradation I

In the field of microgrids with a significant integration of Renewable Energy Sources, the efficient and practical power storage systems requirement is causing DC microgrids to gain increasing attention. However, uncertainties in power generation and load consumption along with the fluctuations of electricity prices require the design of a reliable control architecture and a robust energy management system for enhancing the power quality and its sustainability, while minimizing the associated costs. This paper presents a mixed approach illustrating both simulation and experimental results of a grid-connected DC microgrid which includes a photovoltaic power source and a battery storage system. Special emphasis is placed on the minimization of the total operating cost of the microgrid while considering the battery degradation cost and the electricity tariff. Thereby, an optimal energy management system is proposed for Energy Storage Systems scheduling and enabling the minimization of the electricity bill based on simple models. Simultaneously, the differences between simulation and laboratory performances are highlighted.

INTRODUCTION

With a considerable increase of Renewable Energy Sources (RESs) in electrical energy systems, the concept of Microgrids (MGs) has emerged a novel paradigm for active distribution networks contributing to the emergence of smart grids. The latter concept has been gaining increasing interest by the research community during the last 20 years [START_REF] Vasilakis | The Evolution of Research in Microgrids Control[END_REF]. MGs can be connected to the main electricity distribution grid in order to improve reliability of power production and reduce power losses due to the proximity between energy productions and consumptions as shown in Fig. 1. Nevertheless, they can also be operated in an islanded mode to improve the resilience, ensuring the power supply continuity in case of grid failure as they are potential solutions for electrification of off-grid islands, rural villages, and remote areas [START_REF] Nasir | Solar PV-Based Scalable DC Microgrid for Rural Electrification in Developing Regions[END_REF], [START_REF] Nduhuura | Impacts of Electricity Outages in Urban Households in Developing Countries: A Case of Accra, Ghana[END_REF]. Furthermore, in case of default of the utility grid, some MGs can collaborate with other MGs in the neighborhood. Nowadays, MGs can be fine-tuned so as to minimize their associated energy bills, while consuming more renewable energy and be better integrated in their environment. Moreover, technical adjustments enable to mitigate the control and coordination problems which can originate from inter-temporal operation of Energy Storage Systems (ESSs) and intermittent nature of RESs issues. Thus, an Energy Management System (EMS) becomes indispensable in order to explicitly handle constraints, energetics and economic criteria thereby enabling optimal scheduling of power resources and ESSs in MGs so as to achieve optimal supply-demand balance [START_REF] Wang | Two-stage energy management for networked microgrids with high renewable penetration[END_REF]. EMS is the tertiary layer of MG control operating hierarchically above the real-time control and the power management levels, which are the primary and secondary layers respectively [START_REF] Yamashita | A review of hierarchical control for building microgrids[END_REF]. In the recent years, a number of research papers related to EMS dedicated to MGs have been published. In [START_REF] Saharia | A review of algorithms for control and optimization for energy management of hybrid renewable energy systems[END_REF], a review of several optimization methods of energy management is presented. A controller for EMS called the Resilient Model Predictive Controller (RMPC) based on the Model Predictive Control (MPC) principles is developed in [START_REF] Prince | Resilience in energy management system: A study case[END_REF] for grid connected residential MG in order to minimize the electricity bill and load shedding. Furthermore, in [START_REF] Prince Agbodjan | Stochastic modelled grid outage effect on home Energy Management[END_REF] a Resilient Open Loop Feedback Controller (ROLFC) for EMS is developed using a multistage stochastic optimization (MSO) approach. In [START_REF] Chakraborty | PV-Microgrid Operational Cost Minimization by Neural Forecasting and Heuristic Optimization[END_REF], a Linear Programming (LP) along with heuristic algorithm is proposed for minimization of operating cost of a grid connected MG. In [START_REF] Almada | A centralized and heuristic approach for energy management of an AC microgrid[END_REF], a heuristic approach for energy management of an alternating current (AC) MG is proposed which results in non-optimal solutions. A multi-objective predictive control strategy is proposed in [START_REF] Shivam | A multi-objective predictive energy management strategy for residential grid-connected PVbattery hybrid systems based on machine learning technique[END_REF] for gridconnected photovoltaic (PV) and battery systems for residential buildings aiming to dynamically decrease the cost of electricity compared to traditional systems while reducing the CO2 emissions. In [START_REF] Hossain | Energy management of community microgrids considering degradation cost of battery[END_REF], an economic cost model is used for optimal management of MG Distributed Energy Resources (DERs), taking into account the degradation cost of the cyclic ageing of the battery. However, these EMSs have not been experimentally validated on facilities with actual devices or based on hardware-in-the-loop platforms.

Concerning experimentally verified EMSs, an EMS for battery-based hybrid grid-connected MG is designed and experimentally tested in [START_REF] Luna | Mixed-Integer-Linear-Programming-Based Energy Management System for Hybrid PV-Wind-Battery Microgrids: Modeling, Design, and Experimental Verification[END_REF] which enables the scheduling of optimal power references for DERs and also allows to reduce cost in the MG over a 24-hour ahead forecast data. In [START_REF] Luna | Optimal power scheduling for a grid-connected hybrid PV-wind-battery microgrid system[END_REF], a Mixed-Integer Linear Programming (MILP) mathematical model is proposed in order to optimally schedule the power references of the DERs in a grid-connected hybrid PV-windbattery MG. In [START_REF] Zia | Energy management system for a hybrid PV-Wind-Tidal-Battery-based islanded DC microgrid: Modeling and experimental validation[END_REF], an EMS is developed that optimally balances generation and demand with scheduled power generation of PV, wind, and tidal energy sources with the objectives of maximizing the utilization of RES for an islanded Direct Current (DC) MG. With reference to [START_REF] Zia | Energy management system for a hybrid PV-Wind-Tidal-Battery-based islanded DC microgrid: Modeling and experimental validation[END_REF], the battery ageing cost is considered in detail taking into account several parameters such as temperature, operation and maintenance cost and residual value which the author obtained from the practical regression models of battery degradation developed in [START_REF] Zia | Optimal operational planning of scalable DC microgrid with demand response, islanding, and battery degradation cost considerations[END_REF].

Moreover, a dynamic optimization model is proposed in [START_REF] Moayedi | Unifying Distributed Dynamic Optimization and Control of Islanded DC Microgrids[END_REF] for power sharing among sources in an islanded DC MGs, but the battery degradation cost for sustainable battery lifetime is not taken into consideration. In [START_REF] Parisio | Use of model predictive control for experimental microgrid optimization[END_REF], a MPC approach is applied for achieving economic efficiency in a grid-connected MG operation management. Thereby, three different experiments are executed with different planning horizons in order to assess the feasibility of the MPC-MILP control scheme. A predictive EMS is developed in [START_REF] Michaelson | A Predictive Energy Management System Using Pre-Emptive Load Shedding for Islanded Photovoltaic Microgrids[END_REF] for an islanded PV-battery DC MG for preventing or reducing outage duration. In [START_REF] Xiao | Multilevel Energy Management System for Hybridization of Energy Storages in DC Microgrids[END_REF], a distributed control and non-optimal EMS is proposed for preventing State of Charge (SoC) violations of hybrid ESSs in islanded DC MGs, however optimization model for cost-effective MG operation is not included.

Most of the aforementioned EMSs have focused on economic optimization and have been experimentally validated using hardware-in-the-loop platforms. However, the large majority of researches in the existing literature have ignored the battery degradation cost in the latter models. Nevertheless, relatively few studies [START_REF] Chakraborty | PV-Microgrid Operational Cost Minimization by Neural Forecasting and Heuristic Optimization[END_REF], [START_REF] Hossain | Energy management of community microgrids considering degradation cost of battery[END_REF], [START_REF] Zia | Energy management system for a hybrid PV-Wind-Tidal-Battery-based islanded DC microgrid: Modeling and experimental validation[END_REF], [START_REF] Zia | Optimal operational planning of scalable DC microgrid with demand response, islanding, and battery degradation cost considerations[END_REF], [START_REF] Parisio | Use of model predictive control for experimental microgrid optimization[END_REF] have considered the wear and tear cost of the battery. Moreover, the latter previous works have not explored the difference between simulation and experimental results. As a substantive contribution in the field of power distribution in MGs, this research not only considers the EMS for maximizing economic benefit and optimally scheduling of ESS but it also considers the battery degradation cost and investigates the discrepancy between prototype and real-life performance. Furthermore, this paper explicitly includes some details on the implementation architecture.

In summary, considering all of the aforementioned limitations regarding previous studies on EMS of MG, the following contributions are made in this research:  An EMS strategy is developed, which performs economic dispatch for minimizing electricity cost of a grid-connected Low Voltage Direct Current (LVDC) MG while satisfying the user's demand in terms of uninterrupted power supply. It takes into consideration the degradation cost of the ageing cycle of the battery in order to reflect the true operational cost of a MG.

 The proposed EMS strategy is validated by an experimental set-up. Analyses are carried out to highlight the differences between simulation and experimental performances.

This paper is organized into five sections. Section II illustrates the case study of an LVDC grid-connected MG operation and also details the EMS mathematical modeling. Section III describes the set-up of the MG experimental platform. Section IV presents the simulation and experimental results of the proposed EMS strategy and finally, Section V presents the conclusion of the paper.

II. CASE STUDY

A. System description

A grid-connected residential PV-battery LVDC MG with a dc-coupled architecture is studied in order to enable experimental validation, as illustrated in Fig. 2. It consists of a discrete time model of a PV system which includes a storage capacity for self-consumption by a user who is connected to the conventional grid. The system environment includes two input variables P pv which is the solar production potential and P 𝑙𝑜𝑎𝑑 which is the power demand. The decision variable (P st ) is the power for charging/discharging the battery, the recourse variable (P g ) is the power exchange with the grid and (P batt ) is the internal power of the battery. An EMS is used for supervisory control and a day-ahead optimal power scheduling of the LVDC MG. The Microgrid Central Controller (MGCC) computes the inflow and outflow of the power profile offline, for each flexible entity, namely the storage units and the grid, in order to minimize the cost of electricity over a given day and to ensure power balance in the MG. It is assumed that the model is a deterministic one where P pv and P 𝑙𝑜𝑎𝑑 of the household are already known in advance and also that the conventional grid is always available.

B. Optimization problem

The optimization problem for achieving the technoeconomical goals of the grid connected LVDC MG is mathematically modeled using LP formulation as in [21] instead of MILP algorithms. This choice of a LP in the control allows to face both storage and grid management constraints such as preventing simultaneous charging/discharging of the battery and purchase/selling of energy to/from the grid respectively which are guaranteed cost-wise without using binary variables. The latter variables would have transformed the problem into a MILP formulation that would involve using MILP solvers instead of LP solvers which may result in a longer computing time for larger-scale examples. Hence, in this study LP was preferred in the perspective of future work.

All the variables at the bus level of this generic minimization problem are linked through the power balance relationship at time t which is given as follows:

P g (t) + P st (t) + P pv (t) = P 𝑙𝑜𝑎𝑑 (t)

(1)

The recourse variable P g and the decision variable P st are separated into two components to fulfill the linearization of the system:

P g (t) = P g + (t) + P g -(t) (2)
P st (t) = P st 𝑑𝑖𝑠 (t) + P st 𝑐ℎ (t) (3)

P g + (t)
is the power drawn from the grid whereas P g -(t) is the power sold to the grid. P st 𝑑𝑖𝑠 (t) and P st 𝑐ℎ (t) are the discharging and charging power of the battery respectively.

1) Objective function and Constraints

The objective function that minimizes the cost of the electricity exchanged with the grid is given by [START_REF] Wang | Two-stage energy management for networked microgrids with high renewable penetration[END_REF]. C E is the electricity cost to be minimized over the time horizon t day , which is sampled at the time interval ∆t.

min C E = min ∆t ∑ [c g + (t)P g + (t) + c g -P g -(t) t day t + c st P st 𝑑𝑖𝑠 (t) + c st P st 𝑐ℎ (t)] ( 4
)
Where c g + (t) is the electricity purchase tariff at time t which can vary according to off-or on-peak energy prices. The costs which are constant over the time horizon are: c g -, the selling price of energy to the grid and c st , the battery ageing unit cost associated with charging/discharging decision variables.

In parallel to (1), the following constraints must be satisfied at time t in order to achieve a feasible solution:

P gmax ≤ P g (t) ≤ P gmax + (5)
Where P gmax + ≥ 0 and P gmax ≤ 0 are defined by the contracted powers subscribed by the user.

Consequently, P batt (t) is calculated linearly as follows:

P batt (t) = P st dis (t) η cvs + P st 𝑐ℎ (t).η cvs . η E (6)
Where η cvs is the efficiency of the converter and η E is the energy efficiency of the battery which is more precisely its round-trip efficiency considered at charging mode only. Thus, the boundary of the charging/discharging powers is:

P battchargemax ≤ P batt (t) ≤ P battdischargemax ( 7
)
Where P battdischargemax ≥ 0 and P battchargemax ≤ 0 are related to the battery C-rate which is the measurement of current in which a battery is charged and discharged. In addition, P st 𝑑𝑖𝑠 (t) and P st ch (t) are also bounded:

0 ≤ P st 𝑑𝑖𝑠 (t) ≤ P battdischargemax .η cvs ( 8
) P battchargemax η cvs .η E ≤ P st 𝑐ℎ (t) ≤ 0 (9)
The power exchange with the battery is indirectly constrained through the energy stored in the ESS allowing the calculation of the State of Energy (SoE). The remaining and available energy of a battery is defined by the SoE indicator with reference to its nominal energy capacity and compared to the SoC which estimates what is left in the battery in terms of the amount of charge related to its nominal capacity [START_REF] Mamadou | Definition of a State-of-Energy Indicator (SoE) for Electrochemical Storage Devices: Application for Energetic Availability Forecasting[END_REF]. The SoE is calculated as given below: [START_REF] Almada | A centralized and heuristic approach for energy management of an AC microgrid[END_REF] Where E nom is the nominal energy capacity of the battery. However, the battery has a limited storage capacity:

SoE(t+1) = SoE(t) - P st dis (t)×∆t η cvs ×E nom - P st ch (t).η cvs .η E ×∆t E nom
SoE min ≤ SoE(t) ≤ SoE max (11)
2) Ageing cost of the battery and grid cost The battery ageing unit cost (c st ) is expressed as the cost of exchanging 1 kWh. Therefore, the total amount of energy that can be exchanged by the battery over its lifetime (E lifebatt ), as derived from [START_REF] Hernández-Torres | Techno-economical optimization of wind power production including lithium and/or hydrogen sizing in the context of the day ahead market in island grids[END_REF], is calculated as follows:

E lifebatt = E nom × N cycles × DoD ( 12
)
Where N cycles is the number of cycles and DoD the Depth of Discharge of the battery. Knowing the battery's price (c batt ), 𝑐 st can be calculated by the following equation:

c st = 1 2 c batt E lifebatt (13) A factor of 1 2
is used to account for both charging and discharging modes during a cycle [START_REF] Zia | Optimal operational planning of scalable DC microgrid with demand response, islanding, and battery degradation cost considerations[END_REF]. This is a simple ageing model and the battery price is estimated at 85 €/kWh, leading to a cost of c st equal to 0.235 €/kWh. Moreover, according to the power subscription by the user and based on off-or onpeak energy price in France (regulated tariffs), the actual purchase price of energy c g + (𝑡) from the grid is 0.1360 €/kWh during off-peak hours and 0.1821 €/kWh during on-peak ones. The selling price ( c g

) is constant and equal to 0.10 €/kWh for a PV system ≤ 9 kWc. Nevertheless, as the electricity price has been increasing in the recent years, based on the evolution curve of the electricity tariff [START_REF] Daubaire | Prix de l'énergie : carburants, gaz, électricité, un tiercé en hausse[END_REF], [25], the hypothesis made in this study is that in the future, the purchase tariff for electrical power will be five-fold higher than the actual price and the selling price will increase only two-fold. This extrapolation is clearly hypothetical but the recent cost evolution of energy prices, especially the EPEX-SPOT tariff with selling prices of electricity varying from 0 to 1000 €/MWh makes credible this assumption [25]. Therefore, the purchase electricity costs are considered as 0.68 €/kWh during off-peak and 0.9105 €/kWh during on-peak hours. The selling electricity rate is considered as 0.20 €/kWh. It is to be noted that off-peak hours can be consecutive or scattered within several time slots. There are 8 off-peak hours per day which are concentrated at different times of the day when demand is naturally less sustained. The 8 off-peak hours per day considered in this study is comprised between 12 am and 8 am.

III. SET-UP DESCRIPTION OF THE EXPERIMENTAL MICROGRID

A test bench is developed at the "Laboratoire d'Analyse et d'Architecture des Systèmes" (LAAS) to emulate the operation of a laboratory scale LVDC MG platform as shown in Fig. 3. The latter bench includes a programmable power supply and a loading facility where both consist of integrated converters, thus enabling the emulation of PV production and load consumption respectively. The test bench consists of an ESS which is connected to the DC bus via a bidirectional current converter. Thus, the ESS and the associated converter can operate under real-life conditions involving actual constraints such as variation in the efficiency of the ESS and the converter, change of temperature at the location, etc. The MG is controlled using the DC Bus Signaling (DBS) which is a control strategy where each entity controls the bus voltage level according to the energy exchange in the LVDC MG in order to maintain its power balance [START_REF] Schonbergerschonberger | DC-Bus Signaling: A Distributed Control Strategy for a Hybrid Renewable Nanogrid[END_REF]. Additionally, the above-mentioned bench is driven following a hierarchical control. The primary control is a decentralized structure used to control the converter as has been developed at LAAS [START_REF] Gaetani-Liseo | Impact on Energy Saving of Active Phase Count Control to a DC/DC Converter in a DC Micro Grid[END_REF]. The Local Controller (LC) present on the converter board ensures the above-mentioned primary control. The secondary and tertiary controls are centralized structures incorporated for ensuring the functions of Power Management System (PMS) by the DBS and the EMS respectively. The MGCC communicates via an Inter Integrated Circuit (I 2 C) communication bus with the LC for the converter associated with the ESS. The PV, load and bench instrumentation system communicates through General Purpose Interface Bus (GPIB) communication protocol.

In this study, a pack of lead-acid battery is connected as the ESS of the MG. The battery rated energy capacity is 2.16 kWh and the maximum power for charging and discharging is 0.54 kW at 0.25 C-rate and 0.846 kW at 0.4 C-rate respectively [START_REF]Datasheet YUASA[END_REF]. The energy efficiency of the battery (η E ) is 83% and the converter efficiency (η cvs ) is 96% [START_REF] Gaetani-Liseo | Impact on Energy Saving of Active Phase Count Control to a DC/DC Converter in a DC Micro Grid[END_REF]. The maximum and minimum limits of the battery SoE are defined as 95% and 5% respectively. The maximum power of P g that can be drawn or transferred to the grid is taken as 0.9 kW. The voltage, current and power limit values of the DC bus are given below in Table I.

IV. RESULTS AND DISCUSSION

The PV production and load consumption data for a 24hour ahead period used for the simulation and experimental set-up is extracted from the "Oahu Archive web site" and "IEEE PES Test Feeder" respectively. These latter values are scaled for the purpose of experimental validation. Fig. 4 shows the PV production, load consumption and net power (P net ) which is the power difference between P pv and P load as expressed hereunder: P net = P pv -P load [START_REF] Luna | Optimal power scheduling for a grid-connected hybrid PV-wind-battery microgrid system[END_REF] Fig. 4. PV production, load consumption and net power.

A. Simulation results

The simulation is carried out over a 24-hour prediction horizon. The power references for a LP based optimal scheduling are computed offline for a sampling rate of 10 minutes (Δt = 10 minutes) by the MGCC. The model is developed in Python using the pyomo modeling language [START_REF] Hart | Pyomo -Optimization Modeling in Python[END_REF] and solved using the GUROBI solver. The proposed EMS is tested with an initial battery SoE maintained at 35% as shown in Fig. 5. It is underlined that for the simulation performed the global optimal solution is achieved, as described below. The optimal solution for a minimal cost is achieved by enabling the utilization of the grid during off-peak hours between 12 am and 8 am. The battery is charged between 1:30 pm and 4:50 pm when P net > 0 starting with an initial SoE of 35% to reach a required level of SoE of about 59% so as to be able to meet the user's energy demand during peak hours without drawing power from the grid. It is pointed out that the battery is not charged to its maximum SoE as it is costly for charging and discharging. Also, an SoE of 59% is enough for the time horizon considered. Furthermore, the battery is discharged at the end of the day between 4:50 pm and 12 am so as to avoid drawing power from the grid during low PV production and peak hours. Consequently, the electricity bill is minimized while simultaneously considering the operating cost due to the battery degradation. The electricity cost over the time horizon for this simulation scenario is 0.27€ as per the objective function value. Thus, if the latter value is positive, the user need to pay the electricity bill (deficit) whereas, if the value is negative the user is on profit.

B. Experimental results

The proposed EMS has been experimentally tested at the lab-scale MG of LAAS. Fig. 6 shows the experimental results of the previous scenario for an initial SoE of 35% of the battery. The latter results allow to assess the validity of the offline optimization model used in simulation. Thus, the optimal trajectories obtained offline based on the LP and on the simple models of PV and battery devices are directly injected as the power references in the actual storage system. The PV production and load consumption are also supplied as power references to the test bench for experimental purposes. In Fig. 6, P 𝑛𝑒𝑡 * is the power difference between the PV production and load consumption, P g 𝑒𝑠𝑡 is the estimated power that is drawn and transmitted to the utility grid related to experimental measurement, SoE max = 95% and SoE min = 5% are the maximum and minimum limits of the SoE of the battery respectively and SoE 𝑠𝑖𝑚 is the resultant SoE of the battery with respect to the offline optimization during simulation. Regarding the experimentation process, P batt 𝑊 is the charging/discharging power of the battery measured by the wattmeters enabling the estimation of the battery's SoE est, W . Fig. 7 is a close-up of the experimental results between 1 pm and 12 am. It can be noticed that the experimental results slightly differ from the simulation results, especially the charging/discharging power of the battery. It can be specifically observed that between 1:30 pm and 3:30 pm, the power reference P batt * and the measured power P batt 𝑊 are not the same. It is noted that P batt 𝑊 is unable to follow the power references of P batt * during the charging of the battery. This is because the latter had reached its maximum voltage, thereby requiring to be charged at a constant voltage while simultaneously limiting its current according to the classic Constant Current-Constant Voltage (CC-CV) protocol for battery charging phase [START_REF] Margot | Identification of ESS Degradations Related to Their Uses in Micro-Grids: Application to a Building Lighting Network with VRLA Batteries[END_REF], [START_REF] Reddy | Linden's handbook of batteries[END_REF]. Therefore, experimentally the battery is unable to charge as much as expected during simulation. It is also noted that during the experiment process, in the time range between 11 pm and 12 am, the battery is not fully discharged contrarily to simulation results. This is due to the fact that the actual battery is not charged adequately during the experimental process, contrarily to the charge process forecasted during the simulation. Thus, a proportion of the required energy must be drawn from the grid during on-peak hours to meet the user's demand, causing an increase in the electricity bill. Thereby, for the identical scenario as in simulation, the electricity cost over the time horizon for the experimental set-up is 0.33€. Thus, the percentage error of the operation cost between simulation and experimental is 22.2%. The above-mentioned discrepancy between laboratory experimentation and simulation can be explained by the lack of accuracy of the battery model. In this study, a simple Power-Energy (PE) model is considered for the battery modeling in conformity with most of the EMSs considering an optimization approach as outlined in the literature. This is especially because a linear model derivation must be achieved for LP. Nevertheless, in the existing battery literature, the SoC is the most common indicator used as it gives an accurate estimate of the battery state [START_REF] Reddy | Linden's handbook of batteries[END_REF], [START_REF] Sauer | State of charge -What do we really speak about ?[END_REF]. Thus, SoC is mainly applied in the field of battery, whereas SoE is mostly employed in EMS optimization approach. However, SoC and SoE are significantly different indicators. Thus, in order to increase accuracy of the battery state, SoC 𝑒𝑠𝑡,𝑊 is estimated in the test bench considering the current and the nominal capacity C nom of the battery. On the other hand, the SoE est, W is estimated with respect to the current, the voltage and the nominal energy of the battery.

The SoC gives the state of the battery in terms of capacity expressed in ampere hour (Ah) current-wise. Contrarily, the SoE gives the state of the battery in terms of energy expressed in watt-hour (Wh) power-wise. Thus, as neither the current nor the voltage of the battery are considered in the offline optimization model, the SoC cannot be calculated in simulation. Therefore, SoC 𝑒𝑠𝑡,𝑊 and SoE est, W are estimated [START_REF] Mamadou | Definition of a State-of-Energy Indicator (SoE) for Electrochemical Storage Devices: Application for Energetic Availability Forecasting[END_REF], [START_REF] Reddy | Linden's handbook of batteries[END_REF]- [START_REF] Zheng | Novel methods for estimating lithium-ion battery state of energy and maximum available energy[END_REF] using the values from the wattmeters in the above-mentioned experimental set-up, as follows: [START_REF] Zia | Optimal operational planning of scalable DC microgrid with demand response, islanding, and battery degradation cost considerations[END_REF] The faradic efficiency of the battery (η F ) is 96%; I batt 𝑐ℎ ≤ 0 and I batt 𝑑𝑖𝑠 ≥ 0 are the charging and discharging currents of the battery respectively; V batt is the voltage of the battery. Fig. 8 shows the difference between the SoE of the battery in simulation and experimental results for the initial scenario of an SoE of 35%. The difference between SoE and SoC can also be observed. SoC 𝑚𝑎𝑥 and SoC 𝑚𝑖𝑛 are the maximum and minimum limits of the SoC for the battery respectively. SoC 𝑒𝑠𝑡,𝑊 is the estimated SoC of the battery during the experimentation process. As shown in Fig. 8, the battery is unable to be discharged at the end of the process as it has reached its minimum SoC limit. It can also be observed that the approximate deviation of the SoE between simulation and experimental results is around 3.5% at the end of the test. The latter results and the percentage error of the operation cost confirm the major drawback of the PE battery model. Nevertheless, as per the results of this study, it can be concluded that the errors mainly due to the PE model of the battery are not significant over a 24-hour time horizon. However, the accumulation of the aforementioned errors may have a significant impact for optimization studies based on a longer time horizon. In order to minimize this error, a more accurate battery model can be used which consider both the current and the voltage of the battery in the formulation of the optimization model. Furthermore, for the simulation and the optimization process of the MGCC, the SoE value must be updated in real time with reference to the monitoring system or by using a Battery Management System (BMS).

SoC 𝑒𝑠𝑡,𝑊 (t+1) = SoC 𝑒𝑠𝑡,𝑊 (t) - I batt ch (t).η F ×∆t C nom - I batt dis (t)×∆t C nom (15) SoE est, W (t+1) = SoE est, W (t) - I batt ch (t)×V batt (t) .η F ×∆t E nom - I batt dis (t)×V batt (t)×∆t E nom

V. CONCLUSION

In this study, the total operating cost of a grid-connected DC MG is minimized by an EMS strategy. Emphasis is also put on the cost implication with regards to the battery lifetime. Both simulation and experimental studies are carried out in order to validate the performance of the proposed EMS strategy at a laboratory-scale MG. The results confirm the feasibility of the supervisory control and demonstrate the significant impact of the battery model on the cost function and also influence the SoE. It is highlighted that the percentage error between simulation and experimental results is 22.2% for the operation cost and around 3.5% for the deviation of the SoE. Thus, this study confirms the limitation of the PE battery model. This work can be further extended involving a more precise modeling of the battery. In the next stage of this research, the design of a decentralized controller by a Multi-Agent System (MAS) will be carried out for a more complex case study regarding a Multi-MG network. A realtime optimization will also be implemented and further cases with stochastic description of grid failure will be investigated.

Fig. 1 .

 1 Fig. 1. Schematic diagram of a generic MG.

Fig. 2 .

 2 Fig. 2. Grid connected DC MG architecture.

Fig. 3 .

 3 Fig. 3. Experiment platform of the LVDC MG.

Fig. 5 .

 5 Fig. 5. Simulation results for an SoE of 35%.

Fig. 6 .

 6 Fig. 6. Experimental results for an initial SoE of 35%.

Fig. 7 .

 7 Fig. 7. Close-up of the experimental results for an initial SoE of 35%.

Fig. 8 .

 8 Fig. 8. Comparison of the variations of SoE of 35% between simulation and experimental results relative to SoC limiting values.

TABLE I

 I

	.	GENERAL CHARACTERISTICS OF THE LVDC MG
		Parameter	Values and Units
		[V min bus ; V 𝑚𝑎𝑥 bus ]	[50 ; 60] V
		[I min bus ; I 𝑚𝑎𝑥 bus ]	[-20 ; 20] A
		𝑃 𝑚𝑎𝑥 bus	1 kW