Jean-Paul Caltagirone
email: calta@ipb.fr

Generalized local boundary conditions in terms of accelerations on an edge-based formulation

Keywords: Generalized Boundary Conditions, Galilean relativity, Discrete Mechanics, Conservation of acceleration, Helmholtz-Hodge Decomposition

The objective is to extend the concepts of discrete mechanics based on a local frame of reference to boundary conditions expressed only from the operators of the equation of motion. Taking into account the principle of Galilean relativity extended to uniformly rotating motions, the boundary conditions do not concern directly the velocity but the divergence and the curl of the latter whatever their kind. The gradient of the divergence defined at the extremities of a rectilinear segment of the boundary and the dual curl of the primal curl of the velocity allow to translate locally any kind of boundary condition of the physical domain. Interactions with other segments of the boundary and interior to the domain are defined causally, by the common vertices. The boundary conditions are integrated within the equation of motion in the form of two orthogonal terms of a Helmholtz-Hodge decomposition. The chosen test cases allow to highlight the properties of the boundary conditions that preserve those of the general discrete formulation.

Introduction

Boundary and initial conditions are perceived in mathematics and mechanics as additional conditions to be satisfied to deal with boundary value problems for bounded domains. In continuum mechanics, the equation of motion provides a general solution and the boundary conditions provide a particular solution to the problem at hand. According to the mathematician Jacques Hadamard [START_REF] Maz'ya | un mathématicien universel[END_REF], a well-posed problem is one for which the boundary conditions imposed concern displacements for one part of the boundary and forces for the second. These conditions should lead to certain properties of the solutions, notably uniqueness. Even if the formalism of mechanics has evolved, the conditions always concern the variable, for example the velocity or the displacement and its tangential or normal derivatives at the boundary of the domain. This formalism is closely related to the concept of continuous medium where all variables are reduced at a point and their derivatives expressed by directions within an inertial reference frame.

The mathematical analysis of the compatibility of the boundary conditions and the equation of motion are numerous and perfectly justified since they condition the quality of the numerical simulations performed in solid and fluid mechanics. Most of the recent studies concern the exit conditions of a fluid from a bounded domain satisfying as much as possible the continuity with the external flow generally not known; the comparison with a larger domain allows however to estimate the errors committed and to validate the proposed conditions. The question of the wellposed boundary conditions for the continuous problem has been treated for example in [START_REF] Nordström | The influence of open boundary conditions on the convergence to steady state for the navier-stokes equations[END_REF] in the framework of the Navier-Stokes equations for steady flows. The case of traction exit conditions for unsteady flows [START_REF] Bruneau | New efficient boundary conditions for incompressible Navier-Stokes equations : a well-posedness result[END_REF] brings an efficient solution to solve vortex exit problems. The work [START_REF] Nordström | Well posed boundary conditions for the Navier-Stokes equations[END_REF] [START_REF] Nordström | Well posed problems and boundary conditions in computational fluid dynamics[END_REF] specify the number of boundary conditions that must be associated to the equations in two and three space dimensions in their compressible formulation. The authors propose a procedure to analyze a general time dependent partial differential equation in terms of well-posedness including boundary conditions.

The objective here is not to present a new formulation of boundary conditions for the Navier-Stokes equations but to come back to the notion of continuous medium and its consequences on the way to impose boundary conditions. The principles of equivalence, inertia and relativity stated several centuries ago form the basis of modern classical and relativistic mechanics. The use of these concepts within a global inertial frame of reference, however, attenuates their scope if they are limited to uniform translational motion. According to the classical definitions, the laws of physics must be expressed in the same form in any inertial reference frame, but a uniform rotational motion is considered as accelerated.

The concept of boundary conditions in terms of accelerations differs in spirit from that described by the mathematician Jacques Hadamard for a well-posed problem. The main difference is due to the abandonment of the notion of continuous medium, of those of forces as well as of essential mathematical properties such as the derivation at a point.

The recent notion of discrete mechanics reinterprets the fundamental principles of mechanics, relativity of velocities, equivalence, inertia, ... on the basis of postulates, in particular that of the conservation of acceleration on an elementary segment. The acceleration is thus considered as an absolute quantity measurable at any point and at any time within a local reference frame. The proper acceleration of a particle or a material medium on a segment is equal to the sum of the accelerations applied to it. The global reference of classical mechanics where the two or three dimensional space is associated to coordinates by direction is replaced in discrete mechanics by only two directions of the local reference composed of the two structures, primal and dual. Only two mechanical effects are present, compression and rotation linked to the two longitudinal and transverse celerities. The rigid motions are also two, translation on a rectilinear segment and rotation around this segment. Moreover only two fundamental units are necessary to describe all the quantities of the equations of motion, a length and a time.

In order to make the formulation of the equation of motion perfectly compatible within the domain, the boundary conditions must be adapted in terms of accelerations. These are expressed on the extremities of the segment where the scalar potential is defined and on the segment itself where the rotation acceleration is specified. The selected examples allow to specify the practical implementation of the numerical methodology and to illustrate the quality of the simulations performed.

Physical model

Geometrical structures

The physical model already presented in detail [START_REF] Caltagirone | On primitive formulation in fluid mechanics and fluid-structure interaction with constant piecewise properties in velocity-potentials of acceleration[END_REF], [START_REF] Caltagirone | Application of discrete mechanics model to jump conditions in two-phase flows[END_REF] is exposed here more briefly. The concept of global inertial reference frame is replaced by a local reference frame composed of a primal elementary structure and another dual one on the figure [START_REF] Maz'ya | un mathématicien universel[END_REF]. The physical domain Ω is tessellated by a set of polygons of any number of sides or plane-faced polyhedra to form a primal mesh of n c cells, n f faces, n e edges and n v vertices. The formulation is edge-based, all operators contribute to express the acceleration γ and the velocity v on a segment Γ oriented by the unit vector t. The vectors γ and v exist only on the segment Γ where they are considered constant; they are in fact scalars on an oriented segment of length dh called the discrete horizon. The primitive face S is bounded by the collection of segments Γ defining the positive orientation of it by the unit vector n; the unit vectors are thus orthogonal by construction n • t = 0. Only some meshes have this property, like the one named Marker And Cell [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface[END_REF] of which the methodology exposed here is an extension to unstructured meshes. The dual surface D composed of segments ∆ orthogonal to the primal facets associated with the segment Γ is assumed to be flat and oriented by the vector t. The extremities a or b are the vertices that connect several segments together.

The discrete operators are defined directly from the primal and dual geometric topologies. The gradient ∇φ = (φ b -φ a )/dh is the restriction of the classical gradient to its projection on the segment Γ; the notion of gradient in space has no more sense since φ is only known at the vertices. The primal curl ∇ × v is computed from the components of the velocity v on each segment as a circulation from Stokes' theorem; the velocity in space will remain indeterminate and not necessary. These first two operators are computed exactly from the metrics of the primal structure. The divergence ∇ • v is computed from the fluxes, the velocity on the segment multiplied by the corresponding dual surface and divided by the dual volume; the result is assigned to point a or b. Finally the dual curl ∇ d ψ is obtained using Stokes' theorem by computing the circulation of the vector ψ along the contour of the dual surface D and the result assigned to the segment Γ as a polar vector. These discrete operators verify some properties of the continuous medium in particular ∇ • ∇ d × ψ = 0 and ∇ × ∇φ = 0. This methodology is very close to those developed previously by Shaskov [START_REF] Shashkov | Conservative Finite-Difference Methods on General Grids[END_REF] for the mimetic methods. The physical model is built on the local reference frame defined previously. The conservation of momentum at the basis of continuum mechanics is abandoned in favor of the notion of conservation of acceleration on a segment; this postulate allows to consider that the intrinsic acceleration of an isolated particle or of a material medium is equal to the sum of the accelerations applied to it. Like any vector, the acceleration γ admits a Helmholtz-Hodge decomposition; thus the accelerations imposed can only be of two forms, a gradient of a scalar quantity and a curl of a vector quantity. The principle adopted in this context admits that the law of motion is written in the form:

γ = -∇φ + ∇ d × ψ (1)
The harmonic terms, both divergence-free and curl-free, corresponding to uniform motions are filtered by the discrete operators of this law.

Velocity relativity principle

In all rigor it is not possible to fix a velocity a priori neither theoretically nor for a numerical simulation realized from an equation of motion; otherwise it would mean that it would be an absolute velocity. The two motions of rigid bodies in translation and in rotation must first be strictly filtered by the equation of motion; if this is the case, it is still necessary to transcribe the boundary conditions of the problem posed without having to impose a velocity.

Let us consider the case of a body in circular motion in a geostationary orbit around a mass M . Gravitation exerts on this body an acceleration equal to g = G M/r 2 where G is the gravitational constant and r the distance between the body and the center of mass M . The so-called centrifugal acceleration is exactly equal and opposite and an accelerometer moving with it would detect zero acceleration. In fact it is the Weak Equivalence Principle (WEP) that associates at any moment a gravitational acceleration to the inertial acceleration. Like uniform rectilinear motion, uniform circular motion must lead to invariance of the form of the equations of motion for any combination of these; this is an extension of Galilean invariance and the principle of inertia. The presence of pseudo forces in some equations reflects the fact that they do not respect this generalized invariance.

In the case of a uniform translational motion w 0 added to a motion v 0 such that v = v 0 + w 0 the equation of motion must be invariant which is the case considering that all corresponding terms are filtered by the equation of motion. For a uniform rotation motion such as v = v 0 +ω×r with w 0 = ω × r the divergence of v is equal to the one of v 0 and the primal curl of v is equal to the one of v 0 to a constant 2 ω • n. The inertia terms are written using the inertial potential

φ i = |v| 2 /2 taking into account the equality |v| 2 = |v 0 | 2 + |w 0 | 2 + 2 v 0 • w 0 . The inertial terms corresponding to w 0 , ∇(|w 0 | 2 /2) -∇ × (|w 0 | 2 /2 n) cancel globally.
The gradient and the dual curl of the last term 2 v 0 • w 0 corresponds to the coupling accelerations between the two motions including the equivalent of the Coriolis acceleration. It should be noted here that the latter does not correspond to a fictitious action and that it does not result from a change of reference frame. The abandonment of the notion of global reference frame requires replacing it by a local reference frame whose only interactions with the neighborhood are of cause and effect.

If it is not possible to impose v directly on the segment Γ it is always possible to impose an acceleration γ on it. In this case the knowledge of the velocity at a previous equilibrium instant v o will allow to evaluate the current velocity v = v o + γ dt. Thus the knowledge of a velocity field, including at the limits, can only be envisaged if it is known at a time t o of its evolution. The notion of inertial reference frame is not enough to define a current state, this one is completely linked to its previous state. The velocity v o can be smaller or larger than the celerity c of the medium, these are two disjoint notions. Although the formulation developed here is applicable to compressible flows, the examples treated below will be limited to incompressible motions.

Boundary conditions

We define a segment Γ s on ∂Ω, the surface limiting the physical domain, where are located the acceleration γ s and the velocity v s . The rectilinear oriented segment of extremities a and b is of length dh = c l dt named discrete horizon where c l is the celerity of the longitudinal wave and dt is the lapse of time necessary for the wave to travel this distance.

The figure [START_REF] Nordström | The influence of open boundary conditions on the convergence to steady state for the navier-stokes equations[END_REF] illustrates the actions that define the acceleration γ s of the material medium or of a particle along the segment Γ s . They can only take two independent forms, (i) a direct acceleration due to a difference of a scalar potential q between a and b and (ii) an acceleration induced by the circulation of a current r along the closed trajectory ∆. The physical modeling of the acceleration on the contour of a polygon or the surface of a polyhedron bounding the physical domain ∂Ω is similar to that inside the domain Ω. The difference is that the domain is bounded and the direct actions are limited to the walls of the domain where boundary conditions must be set. The proper acceleration γ s of a particle or of the material medium on a segment of the boundary ∂Ω is written as a Helmholtz-Hodge decomposition:

γ s = -∇q + ∇ d × r (2)
The scalar and vector potentials, q and r, are respectively defined by the divergence and primal curl operators of the velocity v s :

q = dt c 2 l ∇ • v s r = dt c 2 t ∇ × v s (3)
where c l and c t are the longitudinal and transverse celerities; in the case of a Newtonian fluid the grouping dt c 2 t becomes ν the kinematic viscosity because transverse effects do not accumulate and r is the instantaneous shear stress.

Figure 2. Edge segment Γ s where the acceleration γ s and velocity v s vectors are located; the edge potential q a or q b is applied to the ends of the segment and the edge vector potential is defined by r.

The potentials q and r are energies per unit mass of compression and shear transferred to the physical system Ω by the boundaries or generated directly within this domain; for example it is possible to generate a source in the domain or impose a local rotation. The knowledge of these quantities during the evolutions of the physical system contributes to decompose the different energies (propagated, diffused, dissipated, ...) and to evaluate the global conservation. The acceleration γ s , like the acceleration γ on any segment of the primal geometrical structure, represents the variation of the total mechanical energy between the ends of the segment Γ [START_REF] Caltagirone | On primitive formulation in fluid mechanics and fluid-structure interaction with constant piecewise properties in velocity-potentials of acceleration[END_REF]. In order to preserve the principle of velocity relativity the velocity v s on the segment is defined by:

v s = v o s + γ s dt (4)
It is thus imperative to know the velocity v o s on the segment to evaluate it at a time t = t o +dt. The instant t o corresponds to a mechanical equilibrium where all quantities are known and strictly satisfy the equation of motion. In practice these can be identically zero fields. The field v s cannot be imposed a priori because it would be an absolute velocity which is contrary to the adopted assumptions. The boundary potentials q and r as well as the relations (2) and ( 4) are the only elements allowing to define locally the set of boundary conditions of a boundary problem.

The differences with continuum mechanics are significant, (i) the notion of global inertial reference frame (x, y, z) is abandoned and the boundary conditions are no longer expressed with partial derivatives of the velocity per direction, (ii) the imposition of the boundary conditions is realized locally per segment independently of the spatial dimension considered and (iii) they can of course depend on time.

In summary, for a fluid, the discrete boundary conditions are:

           q = -dt c 2 l ∇ • v s r = -ν ∇ × v s v = v o s + γ s dt (5)
Only the first two apply on the segment Γ, at the vertices for the scalar potential q and on the segment for r. The last expression is a compatibility condition using the principle of inertia which can substitute for the second one and restrict the localization of the boundary conditions to the trace of the physical domain boundary. It may be more meaningful in some cases to impose the stress rather than the velocity on a wall; for example the simulation of a turbulent channel for a given turbulent Reynolds number Re t leads to impose the viscous stress on the wall equal to this value.

Discrete equation of motion

The actual potentials φ = φ o + dφ and ψ = ψ o + dψ are defined at time t = t o + dt where φ o and ψ o are the so-called retarded potentials known at time t o by similarity to the retarded potentials of electromagnetism [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]. The increments dφ and dψ represent the updates of the compression and shear energies [START_REF] Caltagirone | Application of discrete mechanics model to jump conditions in two-phase flows[END_REF]. These are directly related to the divergence and primal curl operators of the velocity. The equation of the dynamics in terms of accelerations is first derived on the Γ segment and extended to a multidimensional space through the scalar potentials common to several segments, from cause to effect. This concept leads to a formulation independent of the space considered in one, two or three dimensions; in one dimension of space only the effects of compressibility remain. In general case, the equation of motion is then written:

               ∂v ∂t = -∇ Å φ o + |v| 2 2 -dt c 2 l ∇ • v + q ã + ∇ d × Å ψ o + |v| 2 2 n -dt c 2 t ∇ × v + r ã φ = (1 -α l ) φ o -dt c 2 l ∇ • v ψ = (1 -α t ) ψ o -dt c 2 t ∇ × v (6)
with the edge potentials q = dt c 2 l ∇ • v s and r = dt c 2 t ∇ × v s . The inertia terms ∇ |v| 2 /2 and ∇ × |v| 2 /2 n have specific forms in discrete which cannot be expressed in continuum mechanics [START_REF] Caltagirone | On Helmholtz-Hodge decomposition of inertia on a discrete local frame of reference[END_REF]. Finally, the Helmholtz-Hodge decomposition addresses all the terms which translate the local and global mechanical equilibrium.

The last two expressions of the system (6) are explicit updates of the potentials computed from the divergence and the curl of the velocity. In the incremental time process the quantity φ o -dt c 2 l ∇ • v is replaced by the update φ of this same quantity in the next step and similarly for ψ. The quantities α l and α t are attenuation factors for longitudinal and transverse waves. Like the velocities c l and c t these properties must be simply known and can of course depend on potentials and time. In fluids the propagation of transverse waves is attenuated very rapidly and it is legitimate to adopt an attenuation factor α t equal to unity. In this case the instantaneous shear stress is equal to ψ o = -ν ∇ • v where ν the kinematic viscosity replaces the dt c 2 t grouping. For a gas such as air, which is compressible a priori, the ratio of the kinematic viscosity ν and the velocity c 2 is of the order of magnitude τ s ≈ 1.65 10 -5 /340 2 ≈ 1.4 10 -10 s; for values greater than this time the importance of transverse acoustic waves can be neglected.

It should be noted that in discrete mechanics the scalar and vector potentials are not variables as such, they are energies accumulated over time from an instant of mechanical equilibrium. They are completely subject to the velocity or rather to its divergence and its primal curl. The pressure which is extracted from them in the form p = ρ γ φ, where ρ is the density and γ the ratio of the specific heats, is not a variable either; even if it can be calculated at a posteriori this quantity is not necessary, it introduces at this stage constitutive laws unnecessarily. All the quantities of the system (6) are expressed with only two units, a length and a time; as the mass is a representation of the energy it is redundant to introduce the two notions.

When the motion is steady the terms ∇φ and ∇ d × ψ are equal; as these two quantities of different variables are orthogonal they can be separated by applying to the equation ( 6) successively the operators divergence and primal curl considering that ∇ • ψ = 0:

           -∇ 2 Å φ o + |v| 2 2 -dt c 2 l ∇ • v + q ã = 0 ∇ 2 Å ψ o + |v| 2 2 n -dt c 2 t ∇ × v + r ã = 0 (7)
The integration of one of these forms is only possible if the boundary conditions are known; in the general case q and r are not independent for any flow, what is known is the sum γ s = -∇q + ∇ × r linked to the velocity v s .

The discrete formulation (6) thus corresponds to an entanglement of direct and induced effects as imagined by J.C. Maxwell [START_REF] Maxwell | A dynamical theory of the electromagnetic field[END_REF] by associating the electric and magnetic effects through a derivative in time. Physically, there can only be an exchange between the compression and shear energies (rotation) if the phenomenon is unsteady. This system of equations has obvious symmetries to be related to the properties attributed by E. Noether to the equations of physics composed of Lagrangians or Halmiltonians. The two terms on the right of the equation ( 6) are two Lagrangians associated with compression and rotation effects respectively. The conservation of the total mechanical energy per unit mass is closely associated with the conservation of the acceleration on a segment [START_REF] Caltagirone | On a reformulation of Navier-Stokes equations based on Helmholtz-Hodge decomposition[END_REF]. It is therefore no longer necessary to ensure the conservation of mass explicitly by an associated equation as in continuum mechanics, there is equivalence between the conservation of energy and mass according to the theory of special relativity. The discrete equation ( 6) is both a conservation of the total mechanical energy, a representation of the momentum theorem and of the angular momentum theorem of classical mechanics; indeed the solenoidal component of the acceleration γ is equal to the dual curl of the vector potential,

γ ψ = ∇ d × ψ.
The numerical methodology implemented to perform simulations from the system (6) has been described in detail [START_REF] Caltagirone | Application of discrete mechanics model to jump conditions in two-phase flows[END_REF]. It is very close to the mimetic methods [START_REF] Hyman | The orthogonal decomposition theorems for mimetic finite difference methods[END_REF] [9] [START_REF] Perot | Discrete conservation properties of unstructured mesh schemes[END_REF] and the Discrete Exterior Calculus (DEC) methods [START_REF] Desbrun | Discrete exterior calculus[END_REF] [17] [START_REF] Mohamed | Discrete exterior calculus discretization of incompressible Navier-Sstokes equations over surface simplicial meshes[END_REF]. However, it differs from the latter by using an edge-based formulation that reflects the physical principle adopted, that the intrinsic acceleration of the material medium on this segment is equal to the sum of the compression and rotation accelerations.

Examples

Examples provided below are intended to implement the prescribed boundary conditions in this discrete context. The test cases chosen correspond to situations for which there is an analytical solution or numerical simulation results reported by many authors. They are limited to simple verification or validation tests.

Uniform motions, source and rigid rotational flow

The uniform translational motion is defined by the vector field v = u 0 e x + v 0 e y where u 0 and v 0 are constants. The uniform rotational motion is written as v = ω 0 × r where ω 0 is the rotational velocity assumed constant and r is the distance to the axis of rotation. These two motions are filtered by the operators, indeed the first one is both zero divergence and zero curl and the second one is zero divergence and constant curl and disappears in the sequence of the operator ∇ d × ∇ × v. The inertia terms ∇ |v| 2 /2 and ∇ × |v| 2 /2 n are also zero for these two movements [START_REF] Caltagirone | On Helmholtz-Hodge decomposition of inertia on a discrete local frame of reference[END_REF]. In R 3 these two motions are eliminated from the discrete equation of motion allowing to extend the Galilean principle of relativity of velocities and that of inertia. In the case where the physical domain is bounded by its limit Γ s it is more complex. The operators ∇(∇ • v) and ∇ d × (∇ × v) on the edge of the domain are no longer identically zero. For example, for a uniform translation motion v = u 0 e x , the divergence of the velocity on the edge, in x = 0, represents the flow q which enters the domain and it is then dt c 2 l ∇ • v + q which must be zero. To verify these assertions, two emblematic examples are used, a point source placed at the center of the domain and the rigid rotational flow v = ω 0 × r.

Let us consider the case of a point source defined by the complex potential f (z) = (D/2π)Logz where D is the flow rate of the source and z = x + i y in a square domain [-0.5, 0.5] 2 . The complex velocity ζ = u -i v allows to compute the corresponding velocities in the (x, y) plane. The divergence of the theoretical velocity is not zero but the flow is indeed curl-free. From the information of the flow rate and the velocity on the edge Γ s a simulation was carried out with the system of equations ( 6) with a null initial solution, φ o = 0, ψ o = 0, v = 0. The solution obtained is well at curl-free ∇ × v = 0 and at zero divergence ∇ • (v + v s ) = 0 provided that we consider the velocity v + v s . The figure [START_REF] Bruneau | New efficient boundary conditions for incompressible Navier-Stokes equations : a well-posedness result[END_REF] shows the scalar potential field φ o and the divergent trajectories from the source. Note that the Bernoulli scalar potential φ o B is zero, the pressure increasing as |v| 2 /2 during recompression as the radius increases. The second example takes the rotational motion defined by the field v = ω 0 × r and the function of complex variable f (z) = -(i Γ/2 π) Log(z) where Γ is the constant circulation of the velocity v on any circle (C) around the origin. The divergence of the velocity and the theoretical curl at the edges lead to the boundary values q and r. The solution is sought from the discrete equation of motion integrating this information with zero potential and velocity fields. The resolution obtained after a time step is exact to the machine precision. With ω 0 = 1 the rigid rotational flow is well at zero divergence and at constant rotational ∇ × v = 2 in the whole domain. The figure (3) illustrates the obtained flow where the streamlines are circular and the scalar potential field φ o corresponds to an increasing pressure towards the outside according to the theoretical solution. Contrary to the previous example the solution in this case is exact to machine precision because the velocities have linear variations in space. Thus the uniform rotational motion is strictly with zero divergence and constant curl in the whole domain including for a domain bounded by its limit Γ s .

Poiseuille flow

Steady planar Poiseuille flow

The next case chosen to highlight the application of boundary conditions in terms of accelerations corresponds to the steady incompressible Poiseuille flow in a planar channel shown in Figure [START_REF] Nordström | Well posed boundary conditions for the Navier-Stokes equations[END_REF]. A parabolic Poiseuile velocity profile is imposed on the left wall, the right outlet is maintained at a fixed potential equal to zero and the boundary conditions of the two horizontal surfaces are no-slip. For this flow where the velocity is imposed upstream, the solution is instantaneous because the incompressibility associated with an infinite celerity of the longitudinal waves leads to the establishment of the profile imposed on the whole channel. The initial condition corresponds to zero fields on the velocity v o = 0 and on the potential φ o = 0.

A flux per unit area ∇•v s , with v s = v m (1-y 2 /h 2 )n where v m is the maximum velocity and h the half height of the channel, is imposed on each point of the left boundary. On this boundary the velocity must remain at its zero value which implies that the acceleration on each of its segments must be zero, γ s = 0. We deduce the value of the shear acceleration ∇ d × r = γ s -∇q to be applied implicitly on each of the segments of this boundary in the form of a source term. The local value of r is fixed by that of the velocity profile, i.e. r = dt c 2 t ∇ × v s . The boundary conditions on all the segments of the horizontal walls are defined by a zero acceleration, γ s = 0 and a shear stress which remains that of the input conditions in y = ±h. Thus the potential to be imposed on these walls is not zero, it is equal to φ s such as ∇φ s = γ s + ∇ d × r. The variation of the potential along the horizontal walls is thus linear φ s (x) = -α x and defined to one constant. Finally, the free output condition which is classically written as v s • n = 0 corresponds here to a potential such that γ s = 0, ∇ d r = 0, i.e. ∇φ s = 0 and φ s = Cte. The set of conditions to be imposed on each segment Γ s of the boundary of the physical domain is entirely defined on this trace without any further reference to information inside the domain.

The solution of this problem is obtained in a single iteration with an accuracy close to the machine error from the system of equations ( 6) whatever the spatial approximation used for regular meshes [START_REF] Caltagirone | Application of discrete mechanics model to jump conditions in two-phase flows[END_REF]. The value of the parameter dt c 2 l is chosen so as to obtain a zero divergence at the machine precision. From a technical point of view, the solution of this system of equations is very easy because only source terms are imposed on each segment and its extremities. Contrary to some penalty techniques, the velocity is not imposed very strongly. In a coherent way with the discrete physical model, the solution on the limits requires to know it at the previous time v o s , the principle of relativity of the velocities is thus completely satisfied.

Circular Poiseuille flow leading to vacuum

The inertia terms ∇ |v| 2 /2 and ∇ d × |v| 2 /2 n are not zero for a Poiseuille flow, they compensate instantaneously and locally. The side walls do not exert any action on the fluid and the pressure is constant in each straight section. The mechanical equilibrium must therefore be preserved in the absence of walls which have only a viscous frictional action on the fluid. In the case where the fluid is not viscous, the jet must continue its motion with the same velocity profile according to the principle of inertia.

In order to analyze the mechanisms of this equilibrium in more detail, let us consider the case of a viscous liquid, oil, whose vapor pressure at the temperature considered is assumed to be zero in order to eliminate a possible radial mass transfer. A Poiseuille flow established in a channel with a circular cross-section opens into the vacuum and continues its course following a straight trajectory in the absence of gravity from zero velocity and potential fields. Physically the liquid jet has no reason to destabilize without an imposed disturbance; in the presence of an external gas the oil jet can destabilize and split according to modes established by the linear theory of stability. In the absence of walls, the viscous friction inside the fluid tends to homogenize the streamwise velocity to find the flow velocity and thus ensure the conservation of mass.

The figure [START_REF] Nordström | Well posed problems and boundary conditions in computational fluid dynamics[END_REF] represents the axial and radial velocity fields of the axisymmetric flow where the radius of the upstream tube is R = 0.1 m, the total length is l = 1.1 m, the flow velocity u s = 1 m s -1 ; the kinematic viscosity of the fluid is equal to ν = 5 10 -3 m 2 s -1 . The axial upstream velocity, u ∈ [0, 2], becomes equal to u = 1 and the radial velocity tends towards zero. The inertial effects at the exit of the jet fade away under the effect of viscous constraints to give a laminar steady flow in perfect mechanical equilibrium. The liquid jet remains in its initial geometry and no mass transfer is observed towards the vacuum. Some specificities of discrete mechanics are highlighted during the simulation of this example. First of all, the density of the fluid is not specified because it is not involved in the system of equations that describes the mechanical equilibrium in terms of accelerations; this problem does not depend on the density. However, the formulation used allows the simulation of longitudinal and transverse waves in vacuum for appropriate frequencies. The vacuum as an entity of mechanics corresponds to a field with zero divergence ∇ • v = 0 and zero curl ∇ × v = 0; the acoustic waves are of course represented by longitudinal and transverse accelerations equal to zero c l = c t = 0. In particular a gas whose pressure decreases towards zero becomes infinitely compressible. In the present case the longitudinal velocity is equal to c l = 10 5 in the liquid and to c l = 0 in vacuum, while the radial velocity equal to zero in vacuum becomes the product dt c 2 t = ν in the liquid. The conservation of the mass of the fluid is thus ensured exactly. Taking into account the surface tension per unit mass σ and the capillary potential φ c = σ/R would not change anything on the obtained solution.

For some problems of two-phase flows with free surfaces, the surrounding fluid, generally of lower viscosity than the liquids, does not modify very much the solution obtained by considering it. Its replacement by a vacuum zone can be very effective. However each case is particular and it is necessary to analyze the physics of the phenomena in depth; for example a gas such as air, which is rather compressible, has a celerity c l = 340 m s -1 only 4.4 less than that of water and the role of air bubbles in a water flow is of course not negligible.

Lid-driven cavity

The test-case of lid-driven cavity is of particular interest because the velocity discontinuities at the two upper corners are a source of numerical instabilities which can propagate in a part of the flow. The proposed boundary conditions are thus tested first in the Stokes regime and then for Reynolds numbers of Re = 1000 and Re = 5000.

Stokes flow into lid-driven cavity

Let us first consider the Stokes flow into a square lid-driven cavity of side L = 1 whose walls are of the no-slip type closed by an upper wall animated by a constant velocity v = V 0 e x with V 0 = 1. The upper left corner is represented schematically on the figure [START_REF] Caltagirone | On primitive formulation in fluid mechanics and fluid-structure interaction with constant piecewise properties in velocity-potentials of acceleration[END_REF] where the segment Γ s has for extremities a and b whose velocities are equal respectively to v a = 0 and v b = u 1 = 1; the vertex a is a stagnation point while the velocity on the segment Γ s is, by definition, constant. In fact the velocities on the vertices are only used to define the scalar quantity v 2 a = 0 or v 2 b = 1.

Thus the scalar potential (the pressure) must respect the condition -∇φ i = v 2 av 2 b /dh and be equal to φ o = ∓V 2 0 /2 in the two upper corners of the cavity. Even in the absence of the inertial terms ∇φ i and ∇ d × ψ i the scalar potential (the pressure) is not zero and it is only when V 0 = 0 that this overpressure strictly cancels. In all rigor a Stokes flow only makes sense when v → 0 and in this case the corresponding equation should take the form ∇ 2 v = 0.

Since the velocity initially fixed on the upper wall is equal to v o s = 1 and the acceleration γ 1 is zero it is possible to impose the current velocity v s = 1. On the four walls of the cavity neither the potential q nor the vector potential r is known but the scalar potential q in the corners of the cavity can be obtained by applying the divergence operator to a zero velocity field satisfying only the boundary conditions of the problem. From the values of the potential q thus computed it is possible to obtain the searched behavior in the corners in particular the potential equal to φ o = ∓V 2 0 /2. It should be remembered that the divergence and dual curl operators applied to a bounded domain have a specific behavior at the limits. For example a constant field v = 1 e x has indeed a zero divergence within the domain but the values on the edges depend on the adopted geometrical structure; they represent the local flow of fluid that would have to be injected to obtain a zero divergence on the edge. The figure [START_REF] Caltagirone | Application of discrete mechanics model to jump conditions in two-phase flows[END_REF] represents the scalar potential φ o , the vector potential ψ o and the current function corresponding to this Stokes flow. The solution is free of any parasitic current that could be generated by the presence of discontinuities in the upper corners of the cavity. The perfectly symmetric flow in the absence of inertia presents two Moffatt vortices [START_REF] Moffatt | Viscous and resistive eddies near a sharp corner[END_REF] in the lower corners of the cavity. From the scalar and vector potentials obtained from the numerical simulation it is possible to verify the equality (8) strictly:

-∇φ o + ∇ d × ψ o = 0 (8)
with

ψ o = -ν ∇ × v.
The two vectors of this equality are orthogonal, the gradient of a potential is equal to the dual curl of another potential. We verify by calculating the local primal curl of the velocity that the result allows to find strictly the curl imposed by the upper wall:

Ω ∇ × v dv = Γ v • t dl = -1 (9)
The conservation of the mass ∇ • v = 0 is ensured at the machine precision and the equality [START_REF] Shashkov | Conservative Finite-Difference Methods on General Grids[END_REF] shows the conservation of the curl on the whole domain. The boundary conditions in terms of accelerations are applied to inertial flows corresponding to a lid-driven cavity for Reynolds numbers of Re = 1000 and Re = 5000; these two test cases have been the subject of reference works in particular [START_REF] Botella | Benchmark spectral results on the lid-driven cavity flow[END_REF], [START_REF] Bruneau | The lid-driven cavity problem revisited[END_REF]. The figure [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with a free surface[END_REF] shows the results obtained by two simulations at these Reynolds numbers.

The case of lid-driven cavity at Re = 1000 has been realized using a structured mesh of n 2 c = 16 2 cells and n e = 544 edges tightened to the walls according to a Chebyshev type law. The resolution of the linear systems whose unknowns are the n e edges which support the velocities, is performed using a BiCGStab2 conjugate gradient method or by a direct LU method which is still competitive in dimension two. The flow is clearly under-resolved but the solution is close to that obtained with higher resolutions, in particular we observe the two vortices in the lower part of the cavity. The methodology described for imposing the boundary conditions does not modify the numerical solution, it eliminates the possible perturbations created by the velocity discontinuities in the upper corners of the cavity. The result at Re = 5000 compared to the one obtained previously with the same numerical methodology [START_REF] Caltagirone | On primitive formulation in fluid mechanics and fluid-structure interaction with constant piecewise properties in velocity-potentials of acceleration[END_REF] but with different boundary conditions is almost identical but presents less boundary perturbations on the scalar potential. In order to capture the secondary vortex in the lower right corner it is necessary to increase the resolution of the numerical simulation, it has been increased to n 2 c = 256 2 and n e = 131584 edges. In all the treated cases the divergence of the velocity is of the order of magnitude of the machine accuracy. The use of the Bernoulli potential φ o B = φ o + |v| 2 /2 incompressible flows has some advantages including that of not explicitly calculating the inertia term on the vertices which introduces additional numerical errors; on the other hand, the evaluation of |v| 2 /2 at the center of the faces of the primal structure proves to be much less problematic.

This example is only a verification of the correct implementation of the proposed boundary conditions. If the solutions are the same, as the orders of convergence obtained, the implemented methodology does not introduce additional errors compared to the spatial discretization inside the domain. Above all, it makes the physical modeling of the boundary conditions and the equation of motion in terms of accelerations consistent.

Startup Poiseuille flow in a channel

A Poiseuille flow in a plane channel is initiated from a zero solution v = u e x + v e y = 0 under the effect of a constant pressure gradient δp = 1 e x . The horizontal boundary conditions at y ± e = ±1 are no-slip and the lateral conditions are periodic. The dynamic viscosity of the fluid is equal ρ = 1.

When the time becomes large the solution converges towards the classic parabolic law u(y) = δp e 2 /2 µ 1 -y 2 /e 2 . The transient solution is constructed as the sum of a steady solution and a purely unsteady solution which tends towards zero as time increases. It reads:

             u(y, t) = δp e 2 2 µ Å 1 - y 2 e 2 ã + ∞ n=0 a n (0) exp Å - (2n + 1) 2 π 2 µ 4 ρ e 2 t ã cos (2n + 1) π 2 y e a n (0) = - δp e 2 µ 16 (-1) n (2n + 1) 3 π 3 (10)
The steady solution obtained by solving the discrete equation of motion is exact to machine precision whatever the spatial approximation used; indeed the solution being of the second degree and the spatial scheme of order two makes it possible to obtain an exact numerical solution. The figure [START_REF] Shashkov | Conservative Finite-Difference Methods on General Grids[END_REF] shows the numerical errors in norm L 2 in space and in time. By comparing the numerical and theoretical solutions one notes an order of convergence equal to two in space and time. These results of convergence are verified on any other case of flow governed by the Navier-Stokes equations and the equation of motion of discrete mechanics. Of course errors from other origins, very distorted meshes, unsuitable boundary conditions, too small simulation domain, ... can modify the global convergence of simulations.

Even in the presence of the two non-zero inertia terms in the two-or three-dimensional space equation of motion of the discrete or Navier-Stokes equation, their sum is identically zero and the solution is the same.

Shifted 2D Taylor-Green Vortex

The purpose of this paper is to implement boundary conditions in terms of accelerations when these are arbitrary. The case of the Taylor-Green vortex in dimension two or three in unsteady or steady conditions generally refers to solutions whose boundary conditions are adapted to the physical domain. This is the case of a previous study carried out with the same spatial and temporal discretization [START_REF] Caltagirone | On primitive formulation in fluid mechanics and fluid-structure interaction with constant piecewise properties in velocity-potentials of acceleration[END_REF] where the unsteady solution of a vortex embedded in a cavity shows it is of order two in space and time.

Consider a bounded physical domain (x, y) ∈ [-0.5, 0.5] 2 containing a fluid of kinematic viscosity ν set in motion by a source S adjusted to obtain a non-zero steady synthetic solution from a zero initial condition, v o = 0, φ o = 0. If v 0 is an imposed maximum velocity, the solution in velocity v = u(x, y) e x + v(x, y) e y and in potential φ(x, y, t) is written:

           φ o B = -v 2 0 cos(π(x -x s )) 2 cos(π(y -y s )) 2 u = -v 0 cos(π (x -x s )) sin((π (y -y s )) v = v 0 sin((π (x -x s )) cos(π (y -y s )) (11)
and the corresponding source term equal to:

   S u = -2 V 0 ν π 2 cos(π (x -x s )) sin(π (y -y s )) S v = +2 V 0 ν π 2 sin((π (x -x s )) cos(π (y -y s )) (12)
where x s and y s are the shifts of the solution with respect to the domain coordinates; in the present case these shifts are x s = 0.15, y s = 0.25. In this case the conditions become local for both compression and shear effects. The simulation of this problem in incompressible is realized from a vector equation equivalent to [START_REF] Caltagirone | On primitive formulation in fluid mechanics and fluid-structure interaction with constant piecewise properties in velocity-potentials of acceleration[END_REF]:

∂v ∂t = -∇ φ o B -dt c 2 l ∇ • v + q + ∇ d × (-ν ∇ × v + ψ i + r) + S (13)
where φ o B = φ o + φ i is the Bernoulli potential and ψ i = φ i n. The incompressibility is strongly ensured by choosing a celerity equal to c l = 10 5 which allows to maintain at each step of the computation the divergence at a level close to the machine accuracy. In order to obtain a laminar solution, the viscosity has been chosen equal to ν = 1. The boundary conditions carried by the potentials q and r are obtained from the sought solution [START_REF] Caltagirone | On Helmholtz-Hodge decomposition of inertia on a discrete local frame of reference[END_REF]; the potential q is computed with the divergence of the solution, non-zero on the edge of the domain and the vector potential r by the shear stress. It is possible to impose the velocity on each of the segments initially and to ensure that the shear acceleration is indeed zero during the simulation.

When the problem is treated as unsteady, the equation shows that the compression and shear terms are entangled; these terms cannot exchange energy directly because they are orthogonal. In steady state the time derivative cancels as well as the divergence of the velocity. In the present case we observe a dissociation of the inertial terms and the viscous terms which lead to two separate equations:

   -∇ d × ψ i = -∇ (φ o B + q) -∇ d × (ν ∇ × v + r) + S = 0 (14)
However it is not possible to obtain a solution for each of these equations because they remain coupled by the boundary potentials q and r and are nonlinear. The steady solution of the problem (φ o B , ψ o , v) is thus sought by directly using the equation ( 13) starting from a null solution and by assigning the boundary conditions from the initial time. The figure [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] shows the solution, the Bernoulli scalar potential φ o B , and the two components of the velocity v. The coarse resolution used n 2 = 16 2 allows to show that the solution is free of boundary perturbations. In particular we find the symmetries of the theoretical field, the accuracy obtained is the same as if the solution for a domain with conformal boundary conditions would be shifted in the two directions of space by x s and y s . The convergence in space has been realized for steady fields. The figure [START_REF] Caltagirone | On Helmholtz-Hodge decomposition of inertia on a discrete local frame of reference[END_REF] represents the errors in L 2 and L ∞ norms for the scalar potential and the velocity. The chosen case, although very classical to study the convergence of numerical algorithms, can be transformed as here to highlight the properties of the proposed boundary conditions by breaking the symmetries synthetic solution takes the form:

          
S u = -3 V 0 ν π 2 cos(π (x -x s )) sin(π (y -y s )) cos(π (z -z s )) S v = +3 V 0 ν π 2 sin((π (x -x s )) cos(π (y -y s )) cos(π (z -z s )) S w = 0 [START_REF] Desbrun | Discrete exterior calculus[END_REF] The equation ( 13) is then solved by extracting from the theoretical solution the velocity and the boundary potentials (v s , q, r) for a resolution of n 3 c = 128 3 or approximately n e = 6, 3 10 6 unknowns.

The component on x of the velocity v in figure [START_REF] Maxwell | A dynamical theory of the electromagnetic field[END_REF] illustrates the steady solution obtained by applying the boundary conditions according to the presented method. Thus, in two or three dimensions of space, the methodology developed where the local boundary conditions on accelerations allows to approach varied problems in an automatic way. Only the information on the flow rate fixed on the ends of the segment and the localized shear stress on the segment itself, are necessary to treat a fluid flow problem while respecting the inertia principle.

Conclusions

In discrete mechanics, boundary conditions are applied to each segment of the ∂Ω surface bounding the physical Ω domain. They can only take two forms, a divergence applied to the vertices of the segment and a dual curl fixed on this same segment; they respectively translate the flow injected in the domain and the shear imposed on its surface. This methodology can be extended to all directions of space, taking care to impose both forms per segment.

These boundary conditions are expressed in terms of accelerations and, as within the domain, the intrinsic acceleration on a surface segment is equal to the accelerations related to compressible and viscous effects respectively. The proposal preserves the velocity relativity principle by filtering out uniform translational and rotational motions. The discrete equation of motion with its consistent boundary conditions is thus independent of the chosen local reference frame.

The examples chosen concern viscous fluid flows for simple situations of practical interest; they include some of the verification and validation cases classically used in Computational Fluid Dynamics. They all show that the numerical methodology implemented for the boundary conditions does not generate additional errors to those observed inside the domain.

Figure 1 .

 1 Figure 1. Geometric structures of discrete mechanics; the primal structure (in blue) is composed of a collection of segments Γ oriented by a unit vector t bounded by the ends a or b; these segments form a polygonal planar surface S whose barycenter noted c defines the normal n oriented positively according to Maxwell's rule. The dual structure (in red) has a flat polygonal surface D bounded by segments ∆. The unit vectors are orthogonal by construction, n • t = 0.

Figure 3 .

 3 Figure 3. Scalar potential f (z) = (D/2 π) Logz and, on the right, vector potential f (z) = -(i Γ/2 π) Loz(z) obtained from the discrete equation (6).

Figure 4 .

 4 Figure 4. The boundary conditions noted on the figure express that, on all walls the tangential velocity is zero, the normal velocity imposed on the inlet is the parabolic Poiseuille profile, the horizontal walls are no-slip and the right outlet condition is free. The solution obtained is exact to the machine precision.

Figure 5 .

 5 Figure 5. Circular steady Poiseuille flow leading to vacuum with a radius R = 0.1 m, lengh l = 1.1 m, mean velocity u s = 1 m s -1 , kinematic viscosity ν = 5 10 -3 m 2 s -1 , horizontal (upper) and vertical (bottom) components of velocity.

Figure 6 .

 6 Figure 6. Upper left corner of lid-driven cavity; the scalar potential q is represented by circles on vertices of Γ s edges and triangles are the components of velocity v s or dual curl of boundary vector potential ∇ d × r.

Figure 7 .

 7 Figure 7. Stokes flow in lid-driven cavity performed with n 2 = 256 2 regular cells; from left to right: scalar potential φ o , vector potential ψ o and streamlines; the potentials verify the condition -∇φ o + ∇ d × ψ o = 0

3. 3 . 2

 32 Lid-driven cavity at Re = 1000 and Re = 5000

Figure 8 .

 8 Figure 8. Lid-driven cavity at Re = 1000 with n 2 c = 16 2 cells and at = 5000 with n 2 c = 256 2 Chebychev mesh; Bernoulli scalar potential φ o B and streamlines

Figure 9 .

 9 Figure 9. Startup Poiseuille flow at t = 0.1 s; relative error ε in L 2 norm on velocity in function of vertical spatial approximation n with δt = 10 -5 (left) and in inverse of time step 1/δt with n = 4000 (right).

Figure 10 .

 10 Figure 10. Shifted Taylor-Green vortex where x s = 0.15y s = 0.25, from left to right: scalar Bernoulli potential φ o B and velocity components v x and v y with n 2 c = 16 2 cells and n e = 544 unknowns.

Figure 11 .

 11 Figure 11. Taylor-Green Vortex: convergence in space where n 2 is the number of degrees of freedom for L 2 norm (left) and for L ∞ norm (right); ε is relative errors on the scalar potential φ and the velocity v.

through a choice of the theoretical solution which is non-conforming with respect to the natural boundary conditions. In fact they do not introduce additional errors to those of the basic spatial discretization.

Shifted 3D Taylor-Green Vortex

In dimension three the Taylor-Green Vortex case is used to show the convergence properties of the developed methods but also to study the turbulence decay at high Reynolds number and to compute the corresponding energy spectra. Let us quote as an example [START_REF] Tavelli | A staggered space-time discontinuous Galerkin method for the three-dimensional incompressible Navier-Stokes equations on unstructured tetrahedral meshes[END_REF] where this test case is used to show the high degree of accuracy of space-time discontinuous Galerkin method on face-based non-structured staggered dual grid.

The velocity field v = u(x, y, z) e x + v(x, y, z) e y + w(x, y, z) e z of steady 3D Taylor-Green Vortex in the domain [-0.5, 0.5] 3 is shifted in the three directions of space, x s = 0.15, y s = 0.25, z s = 0.35; the chosen solution is written: