{'fraction_non_alphanumeric': 0.026317197111527146, 'fraction_numerical': 0.005063742533654275, 'mean_word_length': 4.64359025960958, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On primitive formulation in fluid mechanics and fluid-structure interaction with constant piecewise properties in velocity-potentials of acceleration

Discrete mechanics makes it possible to formulate any problem of fluid mechanics or fluid-structure interaction in velocity and potentials of acceleration; the equation system consists of a single vector equation and potentials updates. The scalar potential of the acceleration represents the pressure stress and the vector potential is related to the rotational-shear stress. The formulation of the equation of motion can be expressed in the form of a splitting which leads to an exact projection method; the application of the divergence operator to the discrete motion equation exhibits, without any approximation, a Poisson equation with constant coefficients on the scalar potential whatever the variations of the physical properties of the media. The a posteriori calculation of the pressure is made explicitly by introducing at this stage the local density. Two first examples show the interest of the formulation presented on classical solutions of Navier-Stokes equations; similarly as other results obtained with this formulation, the convergence is of order two in space and time for all the quantities, velocity and potentials. This formulation is then applied to a two-phase flow driven by surface tension and partial wettability. The last case corresponds to a fluid-structure interaction problem for which an analytical solution exists.


Introduction
In fluid mechanics the dominant model is based on the Navier-Stokes equations where the principal variable is the velocity vector V and the conservation of mass is ensured through a scalar, the pressure p. In solid mechanics there are several formulations, in constraints or displacements; each one has its interest but the formulation in displacements, the Navier-Lamé equation, is more efficient to consider as it is easier to calculate the constraints if displacements are known than the opposite. The equations of mechanics are not unified, the Navier-Stokes equation requires to be accompanied by a conservation equation of mass to obtain the pressure whereas for the Navier-Lamé equation it is not required. A number of structural difficulties associated with the Navier-Stokes equation make it difficult to solve in a certain number of cases, in particular for two-phase flows.
The computation of the stresses and displacements in a solid in interaction with a moving fluid in complex situations imposes robust and precise numerical methodologies based on a reliable mathematical formulation. More and more work on these issues is using a monolithic approach to fluid-structure interaction. However, the most delicate problem is posed for the Navier-Stokes equations which includes nonlinearities, which is not the case with the Navier-Lamé equation. The numerical treatment of Navier-Stokes equations in complex situations remains valid as is the theoretical order of the stability and convergence of the schemes used.
The formulation proposed here is based on a physical model different from that of Navier-Stokes. It is briefly recalled later but the details of the derivation of the discrete equations are developed in [1]. The formulation based on the discrete mechanics allows access, at the same time, to the velocity V, compression stresses φ and shear ψ; this formulation (V, φ, ψ) can be written as a projection type algorithm by a splitting or be used as it is.
Discrete geometric topologies are composed of a primal topology and a dual topology associated with specific variables where the scalars are located on the points of the primal geometry and the vectors on the edges of the latter. This vision is the generalization of the staggered Marker And Cell method on structured mesh. The polygonal or polyhedral meshes with any number of facets make it possible to pass over the usual 2D / 3D distinction for other methodologies.
The operators of the primal and dual geometric topologies defined here have obvious analogies with those of DEC (Discrete Exterior Calculus) method. This very elegant formalism based on Exterior Calculus takes its sources from Differential Geometry at the beginning of the last century; it has more recently been presented in a complete way, notably by Desbrun [2] and Hirani [3]. Its applications are numerous in fluid mechanics, elasticity, electromagnetism, etc.
The main difficulties encountered in solving the Navier-Stokes equations with variable physical properties are developed first. Then the bases of the discrete formulation (V, φ, ψ) from the concepts of discrete mechanics and objective advantages are highlighted. Examples of conventional problems of single fluid or two-phase flows are given to show that the solutions are in accordance with the results admitted in the literature. A case of fluid-structure interaction finally shows the merits of a unified approach.

Classical methodologies
The resolution of fluid flows is, in most cases, carried out using the Navier-Stokes equation and that of the solid interactions by the Navier-Lamé equation; these equations are not unified but many decoupled or monolithic techniques exist to solve fluid-structure interaction problems. Since the numerical difficulties encountered are mostly related to the Navier-Stokes equation, let's resume its most common form: In the general case of compressible flows with variable properties the application of differential operators, divergence or rotation on the equation (1) generate many difficulties. The introduction of a constitutive law and an energy conservation equation complicates the formulation for a compressible flow.
The case of incompressible flows where the conservation of the mass is replaced by the constraint ∇ · V = 0 removes a number of difficulties but also potentialities, for example to represent dilatable flows with variable density to less to adopt the Boussinesq approximation. Even if the constraint is strictly verified, requires to adopt ∇ · V = 0, the application of the operators to the neighborhoods at the interfaces of the two-phase flows generate artefact which are due to the modeling of the Navier-Stokes equation. If the viscosity remains a variable quantity in time and space, for example by the addition of a turbulent viscosity, the derivation of this quantity produces instabilities which are also due to the formulation of the viscous term of the equation (1).
The variations in density and viscosity lead to hypotheses, approximations and simplifications aimed at formulating the physical problem posed in an approximate manner. In particular, during the last decades, many works have been focused on the resolution of multimaterial problems and two-phase flows. Even if the methods prove to be efficient and precise, they do not palliate the errors induced by the adopted physical model. All these difficulties have one and the same source, the notion of continuum. Is it legitimate to derive quantities that are physical properties of the media? Of course we cannot deny that there are variations of these properties within fluid flows or multimaterial problems but placing at the same level variables and physical properties has important consequences on the numerical methodologies used.
To ensure the constraint ∇·V = 0 two large classes of methodologies exist, primitive variable pressurevelocity formulations (p, V) and formulations of applying the curl operator to the Navier-Stokes equation including the one entitled curl-and-vector-potential formulation (Ω, Ψ).
Among the methods using the primitive variables let us quote the method of artificial compressibility [4] where the pressure is reactualized starting from the divergence of the velocity computed by a step of prediction; it is an explicit method but does not require a condition at the limit on the pressure. Its implicit versions, augmented Lagrangian methods [5], [6], consist in integrating the constraint into the equation of motion as a Lagrangian; these are a little heavy from a numerical point of view but extremely effective. Direct methods are just as effective, they solve both velocity and pressure within a single algebraic system. These latter techniques, however, require robust and parallelizable solvers; if the case of two-dimensional problems can be treated by LU type direct solvers, the three-dimensional problems require iterative solvers, for example of preconditioned conjugate gradients.
The splitting of the Navier-Stokes equations historically derived by Chorin-Temam idea and called fractional-step technique in 1968-1969 [7] to which the projection methods are nowadays related, is very widely used for the resolution of fluid mechanics equations. J.L. Guermond's review [8] provides a mathematical and numerical point of view of these constant or variable density methods [9]. Many other variants analyze these methods from a mathematical point of view, [10], [11], [12], [13], and explicitly show convergence on numerical test cases. Another point of view is developed in [14], [15], [16] and [17] to treat the problem of the update of pressure and that of the incompressibility of the flow which is considered exclusively kinematic. This method applies to highly variable density flows.

Difficulties related to the term of inertia
All projection techniques involve, at one stage or another of the formulation, simplifications more or less justified; even if, from a practical point of view, they are minor and do not affect the representativeness of the solution obtained, they inhibit any chance of an exact projection.
The divergence of the incompressible Navier-Stokes equation given by Gresho ans Sani [18] with constant density is written: where p is the pressure, ν the kinematic viscosity. If a source term f corresponding to gravitational or capillary effects is added to equation (1), the curl-free part of this term disappears thanks to the divergence operator. The permutation of the divergence and gradient operators of the viscous term leads to a simplified form without this term; a part of the physics contained in the Navier-Stokes equation is thus lost. This confirms that the pressure is a simple Lagrangian that satisfies the constraint ∇ · V = 0. What remains anachronistic is that the pressure, a dynamic quantity, is used to satisfy a purely kinematic constraint. The point of view advanced here is that the quantity which makes it possible to satisfy the constraint of incompressibility must be purely kinematic. To show it let's examine the consequences of applying the divergence operator to the Navier-Stokes equation in more detail.
In the Navier-Stokes equations, the inertial term can be written according to the Lamb vector [19], L = ∇ × V × V. If we are interested in the divergence of this vector as Marmanis, [20], Hamman & al. [21], Lindgren [22] or Kozachock [23] have been, it can be shown that ∇ · L = ·∇ × ∇ × V − |∇ × V| 2 . In this expression, the first term is called flexion while the second one is the local enstrophy. A rigid rotation movement is a local motion enstrophy as the rotational motion is at zero divergence but non-zero curl ∇ × V = 2 ω. In this section, the purpose is to examine the consequences of the expression of the divergence of the vector of acceleration γ.
In the framework of continuum mechanics, the inertial term coming from the material derivative writes V · ∇V or −V × ∇ × V + ∇(|V| 2 /2). By considering one of these two forms of the inertial terms, the divergence operator applied to them provides the same expression: where the term I 2 is the second invariant of the velocity gradient, written in Cartesain coordinates: If u and v are the components of the velocity in a surface then the term −2 I 2 can be represented by ∇u ∧ ∇v where the gradient operator is applied in the planar surface (x, y) and where ∧ is the external product carried by the normal to this surface.
Is it possible to define a pseudo-vector I where components are associated to each surface defined t i , t j , the unit vectors on each axis for the three directions of space in Cartesian coordinates. Then I 2 is defined as the scalar product of I by the normal n = t i × t j to each of the planar surfaces of the system. The first three terms of the second member of the divergence of the material derivative (3) are zero when the latter is equal to zero but the term I 2 is not. The reason why this term must be null is linked to vector conformity; I 2 is a scalar composed of the sum of the contributions defined by direction. The pseudo-vector I can thus be represented by the curl of a potential Ψ and the divergence of this term is necessarily zero ∇ · ∇ × Ψ = 0. Although I 2 is a scalar, each of the components of vector I has to be zero so that ∇ · γ be zero too if the flow motion is incompressible ∇ · V = 0. The constraints I = 0 are presented as compatibility conditions to be satisfied by component. They are similar to those used for defining displacements from constraints in solid mechanics. Therefore residues from the Lamb vector will disturb the full projection of the velocity on a zero divergence field. The term I 2 is second order and numerically weak but the principle of an "exact" decomposition is compromised.

Discrete Mechanics framework
The framework of discrete mechanics is composed of a geometric structure and a physical formulation based on intangible principles known since Galileo. The discrete geometric topology can be represented as the extension of the Marker And Cell method [24] to unstructured staggered meshes. Figure 1. Local frame of reference of discrete mechanics; each planar facet S of normal n whose centroid is at M is generated by rectilinear edges Γ of unit vector t, the trihedron (m, n, t) being direct. A shock or contact discontinuity Σ can be located in c between the points of this primitive geometric topology. The scalar potential φ is fixed on the points while the vector potential ψ is oriented according to n; the γ and V components of the acceleration and velocity are constant on the Γ edge. The density ρ and the compressibility coefficient χ are assigned to the point a or b and the kinematic viscosity is constant over the entire surface S.
The polygonal or polyhedral geometric structure is composed of elements such as that represented in the figure (1). Γ edges with unit vectors t are connected by the points a and b of the primal geometric topology to form planar surfaces oriented along the normal n. The vectors, acceleration or velocity, or rather their components are defined on each edge where they are considered constant. The scalar potential φ is defined at the point while the vector potential ψ is carried by the normal n. The direct trihedron (m, n, t) is the local coordinate system which serves as a basis for the representation of all physical quantities. A Σ shock or contact discontinuity can be optionally located in c on the [a, b] edge.
In discrete mechanics the physical properties are not derivable, they are piecewise constant and assigned to different locations according to the modeled physical effect. The density ρ is localized on the points but it is necessary to define a density ρ v constant on the edge which depends on the position of c; this is equal to ρ v = α ρ b + (1 − α) ρ a with α = [a, c]/[a, b]. Likewise, a compressibility coefficient χ is associated with each point of the primal topology. The kinematic viscosity ν is itself defined at the center of the S facet, but it is supposed to be constant on it.
The local frame of reference of figure (1) allows to discard any global description and the use of an absolute or even Galilean reference. Uniform translation or rotation movements must be excluded from physical and mathematical modeling. All interactions between local frames of reference are cause-andeffect through common points or edges that will transmit information through propagation or diffusion from one element to another. The edge of length d = [a, b] can be as small as necessary without ever being reduced to a point and the angles are preserved in a homothetic reduction. The notion of continuous medium is itself abandoned, which results in the loss of certain properties such as the derivation at a point.

Discrete law of motion
The physical modeling at the base of discrete mechanics is based on Galileo's intuitions: i) the Weak Equivalence Principle on the equivalence between gravitational and inertial effects (equality of the inertial and gravitational mass ) and ii) the relativity of velocity (Galilean reference frame). The fundamental principle of dynamics F = m 0 γ where m 0 is the mass at rest and γ the proper acceleration of the body to which the force F is applied becomes m 0 g = m 0 γ when applied to gravitation. The mass can thus be removed from this last expression to make it an equality between accelerations. The postulate of discrete mechanics consists in extending this form to all the other forces per unit of mass by designating by h the sum of the accelerations applied to the body. The equality γ = h interprets the conservation of the accelerations: the proper acceleration of the particle or material medium is equal to the sum of the accelerations which are applied to it. The conservation of accelerations and energy replace those of mass and momentum of the continuum mechanics.
The discrete mechanics has been the subject of some preliminary work [1] in order to unify the mechanics of solids and incompressible or compressible fluids with or without shock waves. It reproduces, with the same formalism, the results deduced from the Navier-Lamé and Navier-Stokes equations, and those from other domains such as phase change transfer, porous media, and so on.
The law of motion in discrete mechanics is written in the form of a Hodge-Helmholtz decomposition of acceleration: Each component of the γ acceleration can decompose a divergence free and an irrotational one, γ φ = ∇φ and γ ψ = ∇ × ψ; they both act without interaction on Γ because they are derived from orthogonal fields. The equation of discrete motion is derived from the conservation equation of acceleration by expressing the deviations of the potentials φ and ψ as a function of the velocity V. These deviators are obtained on the basis of the physical analysis of the storage-destocking processes of compression and shear energies; the first is written as the divergence of velocity and the second as a dual curl velocity. The physical modeling of these terms is developed in [1].
The vectorial equation of the movement and its updates are written: where φ o and ψ o are the potentials of acceleration γ defined at the instant t o ; these potentials are updated at the moment t to pass from one mechanical equilibrium to another over a period of time dt; the −→ symbol defined this upgrade. The quantities c l and c t are the longitudinal and transverse celerities. For fluids the definition of longitudinal velocity is equal to c l = 1/(ρ χ T ) where χ T is the compressibility coefficient; the transverse celerity of a Newtonian fluid must be replaced by the kinematic viscosity, dt c 2 t = µ/ρ. For solids these parameters correspond to compressibility modules, c l = (λ + 2 µ)/ρ, and shear c t = µ/ρ.
The factors α i between zero and unity correspond to the attenuation of the longitudinal and transverse waves. The attenuation factor α t of an elastic solid equal to unity while for a Newtonian fluid this factor is zero, the transverse waves are attenuated in general on very low time constants (10 −12 s) for fluids.
The variable is the component V of the velocity on each of the edges Γ. We observe that uniform velocity movements V = Cte or uniform rotation V = ω r e θ disappear from the outset of the vector equation of the system (6). The velocity V and the displacement U are themselves upgraded from their value at the instant t o .
Upgraded quantities (φ, ψ) should be advected by the velocity field. This advection, independent of the equation (6), can be realized in Eulerian by the resolution of a transport equation or in Lagrangian. Similarly, if density is not constant, for two-phase flows for example, it must be advected by an appropriate technique.
In the context of the choice of a fixed local reference, the proper acceleration γ can be explained by highlighting the terms of inertia resulting from the discrete mechanics.
The form of the inertial terms differs from that of the continuum mechanics, V · ∇V or ∇(|V| 2 /2) − V × ∇ × V and also corresponds to a Hodge-Helmholtz decomposition. These terms may be partially implicited for better time discretization.

Boundary conditions and orthogonality
Special attention has to be be paid concerning the boundary conditions to be applied in order to obtain an orthogonal Hodge-Helmholtz decomposition. This question is widely debated in textbooks about HHD and especially for its applications to the mechanics of fluids where the constraint of incompressibility is treated by projection methods. The need for orthogonality is discussed by Denaro [25] and by Bhatia et al. [26]. Generally, in fluid mechanics, this is the decomposition of the velocity field that is sought in a free divergent part and an irrotational one; this is extended to the physical modeling of the movement of a fluid or a solid which leads to looking a priori to acceleration in the form of a HHD.
Boundary conditions can relate directly to velocity or its derivatives on the S surface of the Ω domain; in the context of the selected geometrical topology, the orthogonal fluxes at the boundaries are expressed only with the divergence of the velocity and the tangential stresses are represented only with the dual curl of the velocity. Boundary conditions are not linked as in a continuous environment and the concept of well-defined problem is not directly applicable.
Alternatively the boundary conditions can relate directly to the potentials and integrated into the discrete formulation in the form: where D = ∇ · V n S/dt is a flux per unit of time orthogonal to the domain and R = ν ∇ × V t a tangential viscous constraint. This equation is free from additional boundary conditions. From the D and R border operators all types of boundary conditions can be imposed; they can also be used to impose a local source or well or a rotation of a subdomain within Ω.
The legitimate question of orthogonality is discussed here by adopting the notations of figure (1) where n is the unit vector orthogonal to the surface S of the physical domain. Consider the two vectors f = ∇φ and g = ∇ × ψ and writing their inner product < f , g > given by orthogonal decomposition theorems [27], [28]: From the discrete point of view the equalities ∇ h × ∇ h φ = 0 and ∇ h · ∇ h × ψ = 0 are strictly satisfied whatever the polygonal or polyhedral geometrical topologies in the whole physical domain; these properties are then verified on facet and each point of primal topology of figure (1); ψ is collinear to n then ψ × n = 0 and ∇ × ψ is orthogonal to n so that ∇ × ψ · n = 0. This result is verified for all facets S delimiting the polyhedron of the discrete domain wich defines Ω. In both cases we have: The HHD decomposition of the acceleration γ is orthogonal regardless of the BC for the proposed geometric description independently of the compressible or non-compressible nature of the flow.

An exact projection method
If the solutions of the equations of the discrete mechanics are the same as those of the Navier-Stokes equations with constant properties, it turns out that other differences appear when the dissipation term is extracted from these equations or if the divergence operator is applied to the acceleration. We have seen in continuum mechanics that the application of the divergence operator (3) generated terms whose existence was due to the notion of continuum.
Consider now the divergence operator applied to the motion equation The difference between continuum mechanics and discrete mechanics relates to the fact that ∇ · γ = ∇ · (dV/dt). In discrete mechanics, the particular derivative writes and the divergence of it reads By introducing Bernoulli's potential φ B = φ + |V| 2 /2 and applying the divergence operator to the equation of motion it comes: In order to express the partial derivative in terms of differences, V n+1 is called the zero divergence velocity at time t n+1 and its prediction ‹ V. The difference of the Bernoulli scalar potential between two times is denoted The original projection algorithm based on discrete mechanics comes in the form: As can be seen φ ′ B is extracted from the resolution of a Poisson equation with constant coefficients. The pressure, a dynamic quantity, uses the density ρ v on the edge Γ and its computation is carried out completely explicitly directly on the mesh which can be 2D or 3D, structured or not. By choosing to fix the pressure at the first point a to an artibitrary value, we are gradually taking pressure on all the points of the primal mesh.
The jumps corresponding to discontinuities of shock or contact can be introduced directly into the vector equation in the form of a gradient or a primal curl. For example, the capillary pressure jump p c can be implied in the form [[p c ]] = ∇(σ κ) where σ is the surface tension per unit mass and κ the curvature on the Σ interface; if ρ v = α ρ b + (1 − α) ρ a and α = [ac]/[ab] the pressure at b will be written: where the phase function ξ assigned to each of the fluids will be equal to one or zero. The position c of the Σ interface on each edge will simply be known; in practice it can be calculated by the selected interface transport technique, Volume Of Fluid, Level-Set, Front-Tracking, [29], [30], [31], etc. The problem of the residual terms that appears with the application of the divergence operator on the motion equation, in particular the viscous terms, no longer exists in discrete mechanics. Indeed, the viscosity is only present within the dual curl reminding that the properties discrete ∇ h × (∇ h φ) = 0 and ∇ h · (∇ h × V) = 0 are strictly checked for any type of elementary topology.
The term "exact" used in this section refers to a projection method directly derived from the divergence of the equation of motion. It turns out that the numerical differences with conventional projection methods are low with constant physical properties but can generate instabilities when density ratios or viscosities are high.

Discrete formulation properties
In continuum mechanics pressure is the scalar potential of acceleration; it is also the case of the scalar potential φ which replaces the quantity p/ρ v including for flows with variable density. The density ρ v on the edge Γ being constant the product (ρ V) is the momentum, also constant in pieces.
The classical rotational-potential vector formulation reveals the rotational vector Ω = ∇ × V and the vector potential of velocity, V = ∇ × Ψ. But in discrete mechanics ψ is neither the curl of V nor the vector potential of V, this is the vector potential of the acceleration. This quantity is not a kinematic operator, it represents a physical reality ψ = ν ∇ × V is the instantaneous spin-shear stress.
The variables (V, φ, ψ) correspond to a real formulation in primitive variables; the velocity is calculated from equation (6) and the potentials are explicitly upgraded according to only a step.
The advantages of the formulation (V, φ, ψ) can be synthesized into a list: • the treatment of compressible or incompressible fluid flows, with or without shock waves, the computation of the stresses and displacements in the solids and the fluid-structure interaction can be accomplished with one and the same system of equations; • velocity, compression and shear stresses are obtained in the same step; • the method used on many cases of fluid or two-phase flows is particularly robust and precise; it is second order in space and time on velocity and potentials if temporal discretization is carried out with Gear type schemes; • it is applicable for non-constant physical properties, density and highly variable viscosities; in particular the fluid-solid interaction is produced without disturbance to the interfaces; • the interfaces are taken into account without any interpolation, in a sharp way; the primal mesh is plated on the underlying materials or fluids; • all quantities are expressed only with two basic units, a length and a time; • the formulation excludes tensors, so there is no vector gradient to calculate; it is of the ready-to-use type, the discrete operators correspond to the spatial discretization; • the integration of jumps of the physical properties but also of the variables is realized directly in the vector equation of the movement in the form of gradient or dual curl; • the declination of this formulation into a splitting provides the opportunity for an exact projection method by applying the divergence operator to acceleration.
These concrete advantages come directly from a revision of the equations of mechanics and the more fundamental one of the notion of continuous medium. If the one-point derivation of the variables is natural, that of the physical properties introduces many difficulties in the formulation of the equations of physics. The discrete medium in no way diminishes their scope and generalization.

Examples
Numerous examples have shown the validity of the system (6) in fluid mechanics. For example, compressible, incompressible, non-isothermal or two-phase flows can be approximated with DM to recover the standard well known results, in particular the classical analytical solutions of Poiseuille or Couette flows. Reference cases, such as the lid-driven cavity, the backward-facing step or the flow around a cylinder, make it possible to show that the DM model converges to order two in space and time for both velocity and pressure. The flows associated with heat and mass transfers including multi-components are reproduced in a similar way. More complex problems of shock waves, like the Sod tube, phase changes, boiling and condensation [32] are treated in a coherent way by integrating discontinuities within the equations of motion. Two-phase flows with capillary effects, surface tension or partial wetting, are particularly well suited to the discrete mechanics model [16].
In the present form, the system (6) is relatively close to the Navier-Lamé equation associated with the study of stresses and displacements in solids. It differs, however, on several points: more particularly, the discrete formulation is established in velocity and the displacement is only an accumulation of V dt, as the velocity is itself the elevation of γ dt. Numerous examples of simple solicitations make it possible to find the solutions of the Navier-Lamé equation. More complex 2D and 3D problems on monolithic fluidstructure interactions have already made it possible to validate the proposed formulation [33], [34]. The vision of a continuous memory medium makes it possible to treat the problems of large deformations and large displacements in a formulation where pressure stress and shear are obtained synchronously without compatibility conditions. Given the original dissociation between compression effects and rotation, the material frame-indifference introduced by Truesdell [35] is satisfied naturally. The complex constitutive laws can be treated without difficulty and only the physical parameters written in the equation of motion must be known.

The Green-Taylor vortex
The Green-Taylor vortex is a synthetic analytical solution of the Navier-Stokes equations that corresponds to a case of incompressible unsteady flow. This case is often used to compute convergence orders in time and space of numerical methods. It is considered here to show that the discrete formulation makes it possible to find the second-order analytical solution in time and space. More particularly, it should enable us to understand the role played by the different terms of the discrete motion equation and also to characterize how they combine to satisfy operator properties identically.
The equation of the unsteady movement reads: The quantities S i = ∂V/∂t and S µ are suitable source terms that lead to a stationary solution. Parameter r makes it possible to maintain, at any instant, the divergence of the velocity below 1/r; this can be considered as zero throughout the simulation. Since the medium is a Newtonian fluid, the accumulation of shear-rotation constraints is zero ψ o = 0 with the instantaneous potential being equal to ψ = −ν ∇ × V. In fact, the dual curl of the vector potential is equal to source term S µ for this form of the equation. For (x, y) ∈ [−0.5, 0.5], the solution is written: with V = u(x, y) e x + v(x, y) e y and ψ = −ν ∇ × V.
The source term calculated from solution (17) leads to a separation between inertial and viscous terms: The potential of Bernoulli φ B = φ + φ i , where φ i the inertial acceleration, is the equivalent of the Bernoulli pressure. We find that ∇×ψ o = S µ . Adding this source term into the motion equation amounts to imposing a vector potential ψ o depending on time. In discrete mechanics, equilibrium is not required by component and only the addition of accelerations on the edge Γ makes it possible to translate the mechanical equilibrium. Here, the solution of the problem (V, φ, ψ) is obtained directly by the resolution of the vector equation (16).
Convergence in time at order two is obtained using a Gear scheme for the unsteady term, the linearization of the terms of inertia |V| 2 /2 in the form (V n · V n+1 )/2 and the potential upgrade by the expression φ n+1 = φ n−1 − 2 ∇ · V n+1 . Spatial convergence has been studied on a structured mesh from the stationary solution but the results are identical for a pre-defined time or in regular unstructured mesh. The number of points n per direction of space is less than or equal to 2048. The figure (2) shows that the convergence is on order two in time and in space in norm L ∞ on scalar potential and velocity. Figure 2. Green-Taylor vortex: convergence in space (n is the number of degrees of freedom) and in time (δt is the time step) for L ∞ norm; for the scalar potential φ in blue and the velocity V in red.
The use of the vector of Lamb L = ∇ × V × V in continuum mechanics leads to the expression of the Navier-Stokes equation by component. However, this vector makes it possible to represent the inertial term in a planar surface (x, y) for a 2D description. It should be noted that the divergence of the Lamb vector L is not zero, even if its curl ∇ × L = 0 is really null, so that the Lamb vector derives from a scalar potential of the considered problem. In three dimensions of space, the Lamb vector is more difficult to interpret. This leads to different turbulence properties in 2D and 3D. Like all simulations performed in incompressible motions, the results of the discrete model are identical to those of the Navier-Stokes equation.

Lid-driven cavity at Re = 5000
The case of the lid-driven cavity, considered at a sufficiently large Reynolds number, is ideal for testing the legitimacy of the inertial term formulation in the equation of motion. From a general point of view, the solenoidal and irrotational parts of inertia are not easily expressed. It seems advisable to present tangible results to the numericians and CFD specialists who are attached to the inevitable Navier-Stokes equations that are widely used. They have shown their relevance in the immensity and variety of the cases modeled and simulated with them over centuries. Despite the different physical models and the changes brought by the equation of discrete mechanics that replaces the Navier-Stokes formulation, i.e. the treatment of pressure that is transformed into Bernoulli pressure, and the writing of the inertia term in the form −∇ × (φ i n) + ∇ (φ i ), the solutions obtained with both approaches are very close to the reference results obtained previously in the literature for the chosen Reynolds number.
The results of [36] for a Reynolds number of Re = 5000, for which the flow is stationary, are reproduced in table (1). The quantitative comparison concerns the amplitude and position of the vortices generated by the detachment of the flow on the walls of the cavity. Very good accuracy is obtained by comparing the DM formulation and the Navier-Stokes equations to reference [36]. Both formulations provide the same physical solution. A convergence study is investigated for this configuration. It shows a spatial convergence rate of 2 for velocity and pressure. The Bernoulli pressure was used to conduct the evolutions in time and the pressure itself was then extracted from it.

Ref.
ψ max x max y max ψ min x min y min  Table 1. Comparison of the results obtained in mechanics of continuous media [36] and those resulting from the DM formulation presented at Re = 5000 for a Chebyshev mesh with 256 2 cells.
Discrete Mechanics does not bring into question the results obtained with the Navier-Stokes equation. As previously observed, the same solutions are obtained at least with a convergence order and an accuracy that are almost identical to computer errors. However, the most amazing thing is that continuous and discrete models differ on many points. One of these fundamental differences is the use of mass conservation that is always associated with the Navier-Stokes equation. On the contrary, the discrete mechanics equation does not use it explicitly and the mass is always strictly conserved for compressible or incompressible movements: the DM formulation behaves as an autonomous equation that does not require any additional constitutive law. In fact, mass conservation is implicitly integrated into the equation of motion through the term dt c 2 l ∇ · V [37]. The stream lines shown in figure (3) are obtained directly from the vector potential ψ = −ν ∇ × V projected on the surface, ψ = ψ · n.
The lid driven cavity example is a classical flow that shows all the interest of the discrete formulation compared to a continuum type approach. First of all, the solutions with constant properties are strictly the same, whatever the problem dealt with, from the analytical solutions of the Navier-Stokes equations to complex flows whose solutions are obtained numerically. The interest resides in the physical understanding that can be exhibited from the scalar and vector potentials that make the equation of discrete mechanics a true extractor of Hodge-Helmholtz components of acceleration.

Droplet acted upon by gravity
To describe the wetting of a realistic surface, we must overcome a range of obstacles associated with defining, characterizing, and modeling the capillary effects in the presence of three separate mediausually a gas, a liquid, and a solid -which meet at the so-called triple line. Flows involving capillary effects depend on two characteristic properties: the surface tension σ ij between the i-th and j-th media, and the local curvature of the free surface Σ of the immiscible media. Do we really need to introduce another physical parameter to describe the concept of wetting? For the past few decades, the proposed solution has been to observe the contact angle θ between two of the media, typically a liquid surface and a solid substrate. This contact angle represents a static measurement that is deduced from an array of theoretical principles and experimental methodologies. Much has been written about contact angles in the literature, but since the perspective developed here diverges significantly from past work, we shall focus on a presentation that aligns with the discrete approach. When the interface is in motion, past approaches have proceeded by introducing a time-dependent contact angle known as the dynamic contact angle; various observation-inspired laws have been formulated to describe the evolution of angle θ over time under certain specific circumstances.
We shall tackle the problem at a global level from the very start -since we ultimately want to know the acceleration γ of the fluid on the edge Γ, it must be possible to phrase any actions on this edge as a sum of contributions. First of all, there must be inertial and viscous accelerations, a gravitational acceleration, and finally a capillary acceleration γ c . The motion of the triple line depends on all of these effects; it doesn't make sense to define a dynamic contact angle that depends exclusively on time a priori. The model proposed by discrete mechanics is based on the following observation: if the interface Σ with the substrate is allowed to move freely, it will naturally attempt to recover the equilibrium position defined by the static contact angle θ. The dynamics themselves are governed by the acceleration γ. The physical behavior of the interface in the neighborhood of the triple line is determined by the need to establish a mechanical equilibrium from the accelerations of each pair of media. The curvature of the interface in this region is in fact the most effective parameter for describing the tendency of the system to return to equilibrium. If the so-called contact curvature κ c [1] is known, then, although the values of the instantaneous curvature κ depend on each of the various effects, they will tend toward κ c at mechanical equilibrium. For example, in the absence of gravity, a droplet of liquid on a non-deformable solid substrate will tend toward having the curvature κ c that corresponds to the static contact angle θ; at equilibrium, the entire interface will have curvature κ c . The relation between the contact curvature κ c and the contact angle θ is established by considering the geometry of the interface. The key advantage of the proposed model is that it avoids introducing a new parameter, since the curvature is already part of the model of γ c . We can simply require the curvature to be equal to κ c on the triple line and allow the dynamics of the system to be governed naturally by the equations of motion themselves.
The spherical-cap shape is modeled directly and the curvature on the primal topology can be calculated by a formula from differential geometry. In the case of a spherical cap, the curvature is exact (κ = 1/R) up to the machine error. Since the contact curvature is κ c = κ, the velocity is zero, but the capillary pressure is equal to p c = γ κ = σ κ/ρ v , as expected. We will consider the case of the unsteady evolution of a droplet on a planar surface later. Any droplet of liquid on a surface that does not satisfy the equilibrium will necessarily generate motion and oscillations due to the inertial effects that lead it toward the equilibrium state with potential φ c = σ κ.
In the presence of gravity, but without hysteresis or inertial effects, the velocity of the droplet on a surface inclined at an angle of ϕ is constant and the acceleration is zero; the equation governing the mechanical equilibrium in this case may be stated as follows: Gravity introduces an asymmetry between the advancing contact angle and the receding contact angle. As a result, the droplet no longer has the shape of a spherical cap or a circular arc. The solution can only be found by solving aquation (6) with the source term corresponding capillary effects. The notions of contact angle and contact curvature are of course linked -bijectively, in the case of a static equilibrium without external forces. From a physical point of view, the contact curvature is a local and intrinsic property of the surface; this is the maximum curvature that the interface can sustain before a flow is generated to find another equilibrium position. For example, in the case of the lotus effect, the mean surface can be viewed as non-wetting, but this property arises from the high curvature of the interface around the hairs of the leaf. Any gravitationally induced effects are compensated by variations in the curvature, which is essentially constant at a sufficient distance away from the contact points. In the general case, it seems impossible to predict the at rest contact angle in the presence of forces such as gravity.
The static and dynamic equilibria themselves are physically complex and it is difficult to give an answer that is both simple and general. The strategy adopted here is to construct a set of equations of motion that can account for any potential capillary effects, then apply this model to perform direct simulations. The term ∇(σ κ ξ) seems to incorporate all of the desired effects. The examples given below help to justify this perspective.
By using this concept of contact curvature, we can guarantee that the local curvature in the neighborhood of the triple line is equal to κ c at the static equilibrium. Without gravity, in the case of perfect wetting, κ c is equal to zero. For a perfectly non-wetting surface, the contact curvature must be smaller or equal to the inverse of the radius of the sphere that would no longer be in contact with the surface. Again, the desired physical phenomena can only be reproduced if the equations of motion themselves are representative. The contact curvature directly describes the equilibrium or disequilibrium in terms of the resultant of the accelerations on the triple line. Phenomena involving hysteresis are not addressed here, but they can be approached with the same formalism by replacing the advancing and receding contact angles with two distinct contact curvatures.
The triple line represented by the black point in figure (4) belongs to both the interface and the solid wall. We can assign a special value of the curvature to this point, the contact curvature κ c . If this value is the same as the rest of the interface, then, like for the simulation with a structured mesh, the velocity is zero and the overpressure in the droplet is therefore equal to p c = γ κ. If κ c = κ, the forces are no longer balanced on the triple line, and the resultant generates motion that spreads to the rest of the droplet by continuity. In conclusion, whether or not the interface aligns with the mesh, the equilibrium of droplet is guaranteed whenever κ c = κ in the absence of gravity, and the contact angle θ can also be used to describe the static equilibrium. If gravity or some other force is present, the contact angle relation is no longer satisfied, but the concept of contact curvature remains applicable at static equilibrium. Figure 5. Equilibrium of a droplet of radius R = 6.25 · 10 −4 under gravity g = -100 m s −2 from an initial state where the droplet is hemispherical and has zero velocity. The maximum pressure variation in the pressure field is ∆p = 188.9 Pa, and the curvature σ takes values in [835,1574]. Conservation of mass is satisfied up to the machine error.
When gravity is present, the static equilibrium of a rest state is more delicate to describe. Exact solutions are much less common; in the general case, a solution can only be given by solving the equations of motion (6) numerically.
Consider a semicircular droplet initially at equilibrium on a surface with a contact angle of θ = 90Ð® At the initial moment, gravity is increased from 0 to g = −100 m s −2 . The weight of the droplet causes it to flatten out and significantly changes its shape. The curvature, initially equal to κ = 1600, becomes variable along the interface, but the maximum curvature gradually tends toward the imposed value of κ c = 1600 over time. Figure (5) shows a rest state very close to the static equilibrium. Gravity exerts a noticeable effect on the pressure field. Unlike the contact angle θ, the concept of contact curvature κ c is extremely robust. It simply describes the action of any effects applied to the triple line. In the description of a discrete medium, any forces are replaced by accelerations, which provides a number of advantages. Furthermore, we don't need to introduce any other information than is already present in the discrete equations of motion. The concept of curvature, which is already defined at any interface to describe the capillary effects, can also be used to describe phenomena involving the triple line.
The system (6) does not require any boundary conditions, neither on a vector variable such as the velocity component V nor a scalar variable. Every constraint is incorporated implicitly [1] in the same way as partial wetting.

A monolithic approach of Fluid-Structure Interaction
Even if the rheology of the medium is more complex, e.g. viscoelastic fluids, non-linear viscosity laws, viscoplastic fluids, time-dependent properties, and so on, we should still be able to represent its behavior under various types of applied stress. In some cases, the shear-rotation stresses may only be partially accumulated. We can describe viscoelastic behavior by weighting the accumulation term of ψ o by an accumulation factor 0 ≤ α t ≤ 1. Fluids with thresholds can also easily be represented by specifying a value of ψ o = ψ c below which the medium behaves like an elastic solid. Many of the difficulties that are typically encountered in rheologies with non-linear viscosities are no longer an issue with this model.
In discrete mechanics, the concepts of viscosity and shear-rotation are exclusively associated with the faces of the primal topology, where the stress may be expressed in the form ν ∇×V in fluids and dtν ∇×V in solids. As an example, let us examine the interaction of an incompressible viscous Newtonian fluid and a neo-Hookean elastic solid. The stress tensor of an incompressible isotropic hyperelastic material is as follows in the neo-Hookean model: where B = F F t is the left Cauchy-Green deformation tensor. In two spatial dimensions, the Cayley-Hamilton theorem can be used to show that the Mooney-Rivlin model of a hyperelastic material is equivalent to the neo-Hookean model. The velocity of the fluid at the upper wall is shown in black, the velocity of the interface Σ is shown in red, and the displacement over time of the solid U at the interface is shown in blue.
We shall study a problem that was published by Sugiyama in 2011 [38]. Consider an elastic band with an applied shear stress generated by the periodic flow of an incompressible Newtonian fluid. The flow is laminar and periodic in x. Given that there are no compression terms, we can solve the problem in one spatial dimension along the y-axis for y ∈ [0, 1]. Suppose that the upper interface follows the periodic motion V (t) = V 0 sin(ω t), where V 0 = 1 and ω = π. The velocity of the lower surface is kept at zero. The solid occupies the lower part of the domain, and the fluid occupies the upper part of the domain; the position of the interface is y = 1/2. The theoretical solution found by Sugiyama was obtained by separating the spatial variable y from the time variable t; a homogeneous solution is found by considering a basis of Fourier functions on the interval y ∈ [0, 1] and exponential functions on the time interval in each of the fluid and solid domains separately. The sequence of Fourier coefficients can be determined from the coupling at the interface by requiring the velocity and the stress to be continuous. We can find a solution V (y, t) directly from the equations of discrete mechanics (6) simply by imposing the relevant conditions at y = 0 and y = 1. The coupling conditions at the interface, namely the continuity of the velocity and the stress, are implicitly guaranteed to hold by the dual curl operator. The notion of a 2D or 3D space does not exist in discrete mechanics. Instead, the operators define the orientations of the normal and tangent directions within a three-dimensional space. Despite this, the hypotheses of this example enable us to solve along a single spatial dimension. The time step is chosen to be δt = 10 −4 to ensure good overall levels of accuracy; by comparing against the theoretical analytic solution, it can be shown that the numerical solution is second-order in space and time (this result is included in textbook [1]). Figure (6) plots the velocity and the displacement of the interface Σ over time. The velocity of the upper wall is also shown. The solution establishes itself very quickly. After just a few periods, the velocity becomes fully periodic. The velocity profiles are shown until t = 10 s. The displacement of the solid over time may be deduced from the relation U = U o + V dt, where dt represents both the differential element and the time increment δt = dt. Note that the displacement is strongly out of phase with the velocity of the interface.
A selection of the velocity profiles in the y-direction are shown in figure (7) once the periodic regime is fully established. The results converge to second order in space and time. Given the absolute accuracy (of the order of 10 −4 s) obtained using a coarse mesh (n = 32), we can conclude that there is no observable error between the theoretical solution and the numerical solution.
One advantage of the fluid-structure interaction for a neo-Hookean model described by Sugiyama is that it has a theoretical solution. This allows us to compare the numerical solutions that we obtain more precisely, but also allows us to develop new concepts, as we did for discrete mechanics in this section. Sugiyama obtained a first-order error in the L 2 and L ∞ norms, whereas the model (6) achieves second-order results with much lower absolute errors. This improvement is ultimately attributable to the separation of the properties at the interface, as well as the fact that no interpolation is performed, despite a fully monolithic and implicit treatment of the fluid-solid coupling. Fluid-structure interactions on 2D or 3D geometries with a moving interface can of course also be solved using the system (6). However, without an analytic solution for comparison, there is little benefit in doing so, since the errors of the various methodologies accumulate over each step of the process. Other more complex constitutive laws can also be modeled. However, although specific research in this domain is interesting in its own right, this does not offer any additional validation of the discrete model. The complete separation of the laws of motion and the constitutive and state laws in discrete mechanics will likely enable us to apply any kind of law to problems such as this one.

Conclusions
The original form (6) of the discrete motion equation corresponds to a physical model different from the one leading to the Navier-Stokes equation. It is based on the Hodge-Helmholtz decomposition of the acceleration into a component with divergence free and another irrotational. The velocity, the variable of the vector equation, makes it possible to upgrade the scalar and vector potentials which physically represent respectively the compressive and rotational-shear stress. This decomposition was possible considering that the physical properties are not differentiable quantities in the sense of the continuum. Considering these properties piecewise constant still allows any problem of fluid mechanics or unrestrained solids to be addressed. Two-phase flows with highly variable properties, the elasticity associated with complex constitutive laws, the fluid-structure interaction are accessible by solving equation (6).
The declination of the formulation (V, φ, ψ) in the form of a splitting has conducted to an original projection method. This projection method derived from the discrete formulation is described as exact insofar as it results from the application of the divergence operator to the equation of motion to give a Poisson equation with constant coefficients (13) without any hypothesis or simplification. It has the same disadvantage as the other projection methods, the existence of a boundary layer near the walls, but it also has several decisive advantages over the other methodologies like for multiphase flows with high density ratios. This technique is order two in time like the formulation (V, φ, ψ).