Application of discrete mechanics model to jump conditions in two-phase flows
AApplication of discrete mechanics model to jump conditions intwo-phase flows
Jean-Paul Caltagirone
Université de BordeauxInstitut de Mécanique et d’IngéniérieDépartement TREFLE, UMR CNRS n o [email protected] Abstract
Discrete mechanics is presented as an alternative to the equations of fluid mechanics, inparticular to the Navier-Stokes equation. The derivation of the discrete equation of motion isbuilt from the intuitions of Galileo, the principles of Galilean equivalence and relativity. Othermore recent concepts such as the equivalence between mass and energy and the Helmholtz-Hodgedecomposition complete the formal framework used to write a fundamental law of motion suchas the conservation of accelerations, the intrinsic acceleration of the material medium, and thesum of the accelerations applied to it. The two scalar and vector potentials of the accelerationresulting from the decomposition into two contributions, to curl-free and to divergence-free,represent the energies per unit of mass of compression and shear.The solutions obtained by the incompressible Navier-Stokes equation and the discrete equationof motion are the same, with constant physical properties. This new formulation of the equationof motion makes it possible to significantly modify the treatment of surface discontinuities, thanksto the intrinsic properties established from the outset for a discrete geometrical descriptiondirectly linked to the decomposition of acceleration. The treatment of the jump conditions ofdensity, viscosity and capillary pressure is explained in order to understand the two-phase flows.The choice of the examples retained, mainly of the exact solutions of the continuous equations,serves to show that the treatment of the conditions of jumps does not affect the precision of themethod of resolution.
Keywords
Discrete Mechanics; Hodge-Helmholtz Decomposition; Navier-Stokes equation; Two-PhaseFlows; Mimetic Methods; Discrete Exterior Calculus _____________________________________________________________________
This article may be downloaded for personal use only. Any other use requires prior permissionof the author and Elsevier Inc.J-P Caltagirone, Application of discrete mechanics model to jump conditions in two-phaseflows, Journal of Computational Physics, 2021, https://doi.org/10.1016/j.jcp.2021.110151. _____________________________________________________________________
The Navier-Stokes equation in its various formulations forms the indisputable base of fluidmechanics, both from a physical point of view and from its ability to provide reliable predictionsthrough its numerical solution. Discrete mechanics does not call into question the validity ofthese equations, but can be presented as an alternative whose relevance can only be evaluated inthe light of a long test. The case of two-phase flows is one of the relatively recent fields coveredby computational fluid dynamics.Two-phase incompressible flows have particular specificities, such as the presence of interfacesformed by several fluids, generating capillary effects or even very large variations in physical a r X i v : . [ phy s i c s . f l u - dyn ] J a n roperties. One of the problems generated by variable interface flows over time is that ofadvection. There are many ways of transporting these surfaces, such as Volume Of Fluid, Level-Set, Front-Tracking, Arbitrary Lagrangian Eulerian, or combinations of these methods. Manyarticles highlight the defects and qualities of each [1], [2]. The problems of solid inclusion flowalso feature in the very abundant literature on theoretical formulations and their implementation.Numerical treatment of discontinuities has naturally been addressed for highly compressible flows,particularly by Fedkiw [3] and many other authors using the Ghost Fluid method [4]. The caseof discontinuities related to scalar equations is also addressed by this same technique [5], [6].The treatment of discontinuities for two-phase flows makes extensive use of this method [7].Other popular numerical approaches consist in using a cut cell method on a Cartesian meshfor incompressible flows [8] or as an alternative to traditional boundary fitted grid methods forcompressible flows with shock-waves [9].The range of difficulties encountered during two-phase-flows simulations is reflected in variouserrors: (i) advection of interfaces, (ii) calculation of curvatures, (iii) location of the interfacein the mesh, (iv) assignment of physical properties by interpolations on the stencil, and (v)implementation of jumps in the numerical formulation. The objective is to present a formulationand numerical treatment of jumps in physical properties, density and viscosity, as well as capillaryeffects directly associated with two-phase flows.The formal framework is that of discrete mechanics [10], which can be presented as analternative physical model of the Navier-Stokes equations. The associated discrete formulationabandons the notion of continuous medium in favor of an equation on bases of differentialgeometry that are slightly comparable to the DEC (Discrete Exterior Calculus) methods [11]or mimetic methods [12]. Each physical effect of the equation is represented by a solenoidalterm and an irrotational one following a Hodge-Helmholtz decomposition. The jumps in thephysical properties and those resulting from the capillary effects are also represented in the sameway. Unlike other techniques, jumps are located within a single mesh on the one hand and arecompletely implicit within the equation of motion. In particular, no additional ghost point isneeded.To present the advantages of the theoretical formulation, the solutions of degree lower orequal to two are used to validate it. As the methodology of discrete mechanics is of order two inspace, it allows precise simulations to be carried out to the machine error, whatever the regularmeshes adopted and the number of degrees of freedom. A second-order accurate method appliedto a complex problem can cover multiple errors, the point here being to verify that the discreteformulation is free of all artifacts . This section defines the formal framework for deriving the equation of motion. Notionssubstantially different from those of continuous media must be specified. First of all, the Galileanor classical inertial frame of reference is abandoned and replaced by a local frame of referencewhere the interactions are of cause and effect, i.e. the information propagates with a celerity c from one local frame of reference to another. It is then no longer possible to change the frameof reference in the usual sense of the term because the celerity is not a constant quantity, itvaries according to the medium. The concept of continuous medium disappears and the physicalquantities, characteristics and variables are no longer defined at a point; they are located on2ifferent positions of the local frame of reference, the choice of which is dictated by the coherenceof the physical effects described.Similarly, a vector is no longer defined by its Cartesian components in a global frame ofreference ( x, y, z ) . If W is a vector of R , discrete mechanics defines the component of W bothas a scalar on a segment Γ directed by the unit vector t and as a vector V = ( W · t ) t . Onlythe quantity V is considered in the formulation; vector W will remain unused even if it can bereconstructed from the components in a three-dimensional space, but also on each of the facetsof the primal geometry. Thus, the gradient operator of a scalar is not the vector of space in theusual sense but its projection on the segment of unit vector t ; it will be the same for the otherdifferential operators used. The notion of gradient of a vector is non-existent, like all tensorsof order equal to or greater than two in classical mechanics. The concept of tensor created forand by mechanics comes directly from the observation of media having characteristics whichdepend on the direction considered, wood or quartz for example. This natural 18th-centuryidea to transpose the directional characteristics of environments to the modeling of mechanicaleffects can be discussed and modified. In fact, the adoption of a ( x, y, z ) three-dimensional spacedescription three centuries ago persists in present-day mechanics and physics, for example in thetheory of relativity with a four-vector concept. The abandonment of the concept of continuousmedium entails that of derivation, integration and analysis. Discrete mechanics reconstructsa discrete physical model based on simple operators applied to scalars, the vertices and thebarycenters of the facets associated with the normals of the primal geometry. Figure (1) specifiesthe elementary discrete geometry, constituted by a rectilinear segment Γ of ends a and b andof planar facets formed by a collection Γ ∗ of segments delimiting polygonal facets S having acommon side with the segment Γ . Figure 1.
Discrete geometric topology: a set of primitive planar facets S are associated with thesegment Γ of unit vector t whose ends a and b are distant by a length d . Each facet is definedby an contour Γ ∗ , a collection of 3 segments Γ , is oriented according to the normal n such that n · t = 0 ; the dual surface ∆ connecting the centroids of the cells is also flat. The unit vector n located at the barycenter of each facet is orthogonal to that carried by thesegment Γ , the vector t , n · t = 0 . The surface ∆ is not, in the general case, necessarily planarif the polygons delimiting the facets are not regular. This point will be discussed during theprocessing of the jump conditions but, for this step establishing the physical model, this will bethe case. Thus the facets will be regular polygonal, equilateral triangles, quadrangles, hexagons,etc. The physical domain will be tessellated by the pattern represented by figure (1), the primalmesh. 3 primary physical meaning must be attributed to the elementary stencil of figure (1); itis more easily understood in elementary electromagnetism: an electric current in a rectilinearconductor produces a magnetic field in its vicinity and conversely a magnetic field produces anelectric current. Thus a flow represented by a velocity on Γ produces a field of rotation, the curlof velocity ∇ × V , an axial vector carried by the unit vector n . The variation of the velocityon the segment Γ is thus produced by two fundamental actions: a direct action generated by apotential difference between a and b and another action induced by the circulation of the axialvector on the contour of ∆ . These two contributions of the acceleration have no interactionsbetween them because they are associated with orthogonal quantities. They can only exchangeinformation, energy for example, if the phenomenon is dynamic; J.C. Maxwell [13] expressed thisdynamic duality very well in the pathway which led him to the equations which bear his name.The acceleration denoted γ is likewise a component of the acceleration vector, knowledge ofwhich is not useful for physical modeling or the derivation of the equation of motion; like velocity,it will be constant on the segment Γ . In discrete mechanics it has the status of absolute quantityand can be measured without any reference to the external environment; the other quantities -velocity, pressure, energy, potentials, etc. - are relative and defined up to constants which mustbe filtered by the discrete operators in order to ensure the invariance of the equation of motion.Space and time are linked by the local celerity of the propagation phenomenon considered(swell, sound, light) through the notion of discrete horizon dh , such that dh = c dt where dt isthe elapsed time between two states of mechanical equilibrium. The discrete horizon is of theorder of magnitude of the length d = [ a, b ] of the segment Γ . The state of mechanical equilibriumis itself defined as the exact satisfaction of the equation of motion; t o is the reference state and t = t o + dt is the current state. The equation of motion must allow the calculation of the set ofvariables at time t from that at time t o in an incremental time process. The physical model is developed from the following principles and postulates:• Galileo’s remarkable intuition on the principle of weak equivalence between gravitationaland inertial effects of mass. The equivalence between “the grave mass” and “the inertialmass” [14] is now verified at less than one part on , hence its status as a principle;• the Galilean principle of velocity relativity, which should lead to an invariant physicalmodel for uniform motion;• the principle of equivalence of energy and mass formulated by A. Einstein, resulting fromwork on the theory of relativity;• the Helmholtz-Hodge decomposition where any vector can decompose into a divergence-freecomponent and another curl-free component.The principle of equivalence stems from the experiments of the beginning of the 16th century,in particular that of Galileo where two different masses fall with the same acceleration andthe same velocity; Galileo himself attributes this observation to the fact that the mass relatedto inertia is equal to the mass attached to gravity. The weak equivalence principle (WEP)states that the inertial mass and the gravitational mass are equal. Whatever the forms adoptedthereafter, i.e. the equivalence principle of general relativity, strong equivalence principle etc.,presented as a local equivalence between gravitation and acceleration, mass is always presentthere. Singularly, while the phenomenon does not depend on mass, this concept still persiststoday where mass is most often attached to acceleration or velocity. Discrete mechanics returnsto the interpretation of this principle based on the original observation, with the equality of the4ravitational acceleration and the intrinsic inertial acceleration thus expressing the instantaneousmechanical equilibrium of a body.The fundamental law of dynamics or Newton’s second law translates, in its modern version,the equality between the variation of the momentum of a body and the sum of the externalforces, i.e. m γ = F . In this presentation, as the velocities are much lower than the celerity oflight, the moving mass m will be assimilated to the rest mass m . When the force is linked togravity, the law of dynamics becomes m γ = m g or γ = g , an equality between accelerations.The laws of mechanics are still based today on the conservation of momentum q = m V ; influid mechanics the Navier-Stokes equation is a law of conservation of momentum. In physicsthe theory of relativity also considers momentum in the quadrivector formulation. The upkeepof mass or density in the equation of motion was due, until the 20th century, to the legitimatelypredominant role of gravity. What is true for gravity could be extended to other accelerations.In discrete mechanics the fundamental law of dynamics translates that the intrinsic accelerationof a material medium on a segment is equal to the sum of the accelerations which are applied toit: γ = h (1)where h is the sum of the accelerations, those of the effects of compression and shear but alsoall the other potential source terms: gravitation, capillary acceleration, etc.It is essential for the equation of motion to be invariant with respect to uniform motions and,of course, the uniform translational motion at constant velocity; this is Galilean relativity. Theinvariance must be extended to the uniform rotational motion at constant velocity Ω such that V = Ω × r . In discrete mechanics the rotational invariance is structurally assured [15]. Noether’stheorem [16] establishes the link between invariances and the conservation of certain quantities.The rotational invariance leads to the conservation of angular momentum.The equivalence between mass and energy, one of the consequences of the special theoryof relativity carried by the law e = m c , induces a simplistic interpretation which consists inasserting that the mass is equal to energy; however, the velocity v of a particle or a materialmedium is not always equal to the celerity of the wave, whatever its nature (swell, sound, light).The energy per unit of mass is neither equal to v nor equal to c , the equation of motion willgive its expression. However, it is certain that the two concepts are concurrent and it is possibleto conserve mass or to conserve energy, but the conservation of both is redundant. In all theunits of physics, mass only intervenes in the order , or − and it is therefore possible to defineall these quantities per unit of mass. Discrete mechanics establishes that the equation of motionis both a law of conservation of acceleration and a law of conservation of energy.Finally, the Helmholtz-Hodge decomposition commonly used in mathematics, for exampleto project a field onto a space with zero divergence for the [17] projection methods, has onlya marginal role in physics. The Helmholtz-Hodge decomposition makes it possible to write avector in a divergence-free component and another curl-free one, but decomposing the velocityvector in physics involves the relativity of this quantity, which results in the presence of a thirdharmonic component, both divergence-free and curl-free. The Helmholtz-Hodge decompositionof acceleration is of an entirely different nature because the uniform translational and rotationalmotions are filtered out by the discrete equation of motion. The law of dynamics thus becomes: γ = −∇ φ + ∇ × ψ (2)where φ is the scalar potential and ψ is the vector potential, these are the current potentials (atthe instant t ) of the acceleration; they are expressed in m s − , the unit of energy per unit ofmass. The first term of the right-hand side is an acceleration produced by a direct action dueto a difference in scalar potential between a and b , while the second is an acceleration due to an5nduced action generated by the circulation of the vector ψ along the contour of the dual surface ∆ . The derivation of the equation of motion is performed on the segment Γ where the componentsof the three terms of the law (2) are expressed where they are constant. (cid:90) Γ γ · t dl = − (cid:90) Γ ∇ φ · t dl + (cid:90) Γ ∇ × ψ · t dl (3)The length d of the segment Γ is chosen according to the problem considered and is as smallas necessary, but can in no case be reduced to zero; its orientation will thus be preserved. Theintegration of an acceleration on a segment corresponds to an energy, and that of the intrinsicacceleration is the sum of the compression and rotation energies: Φ b − Φ a = (cid:90) ba γ · t dl (4)where Φ is total energy per unit mass. This section is dedicated to the interpretation of discrete vector fields, considering the lineintegrals along the oriented segments of the geometric topology of figure (1). The objective is todeduce a local shape of the continuous medium type. The notations are particularly confusingfor the primal and dual curls and do not allow one to distinguish primal curl ∇ × W , an exactoperation, from dual curl ∇ × ψ , an approximate operation. The geometrical reduction towardsa local equation is carried out differently, however, by considering a homothetic transformationpreserving the angles, in particular. However, this continuous version has the advantage ofallowing a direct comparison with the equation of motion of continuum mechanics. The discreteversion will be established in section 2.2.The discrete equation of motion is thus presented as a law of conservation of acceleration butalso of mechanical energy; acceleration γ , a component of the acceleration vector, is the materialderivative of the velocity component on Γ : d V dt = −∇ ( φ o + dφ ) + ∇ × ( ψ o + d ψ ) (5)The law (2) shows the current potentials while at instant t o only the potentials φ o and ψ o areknown; they are called retarded potentials, like those of Liénard-Wichert in electromagnetism[18]. During the time period dt , the potentials vary in quantities dφ and d ψ and must be modeledfrom the longitudinal and transverse celerities c l and c t , which are different for an elastic mediumbut equal c l = c t = c for a fluid, where it is commonly called the speed of sound. The physicalmodeling of dφ and d ψ are described previously [19, 10].The modeling of the variations of the potentials, dφ and d ψ , can be explained simply byconsidering that these quantities are energies per unit of mass which evolve during a mechanicaltransformation. The first corresponds to an elementary longitudinal compression of a medium ofcelerity c l , where the variation of energy is written dφ = c l ∇ · U where U is the displacement inthe direction considered. The second variation, that of the shear energy, orthogonal to the first,is d ψ = c t ∇ × U ; the latter is polarizable in a plane orthogonal to direction n .Intrinsic acceleration γ is then the sum of two contributions where the first, known as direct,is due to the difference in scalar potential between the two ends of segment Γ . The second isan induced effect of the first; indeed, the velocities carried by contour Γ ∗ create a curl carriedby normal n to facet S ; this axial vector ψ , associated with that of all the other facets having6egment Γ in common, makes it possible to calculate a circulation on the contour of the dualsurface ∆ passing through the barycenters of the facets. This dual curl is carried by segment Γ ,thus defining a second contribution to the intrinsic acceleration γ : (cid:90) Γ γ · t dl = − (cid:90) Γ ∇ (cid:0) φ o − dt c l ∇ · V (cid:1) · t dl + (cid:90) Γ ∇ × (cid:0) ψ o − dt c t ∇ × V (cid:1) · t dl (6)The second contribution is orthogonal to the first and they cannot combine, but simplyoverlap; exchanges between the two contributions are only possible if acceleration γ is not zero.This phenomenon is in line with J.C. Maxwell’s idea on the dynamic role of electromagneticinteractions [13]. Equation (5) is a formal Helmholtz-Hodge decomposition of acceleration.The third harmonic term at the same time with divergence-free and curl-free does not existhere because it corresponds to the uniform translational and rotational movements immediatelyeliminated by the operators in accordance with the principle of relativity. The intrinsic acceleration γ = d V /dt has a particular status: it is the only quantity that can be considered as absolute. Aswe perceive that the direct actions represented by ∇ φ and induced actions given by ∇ × ψ areentangled, in several dimensions of space one does not exist without the other – one is the dualof the other. The classical interactions of electromagnetism between direct currents and inducedcurrents are of exactly the same nature in mechanics. The variation due to an acceleration ofone of the effects implies the variation of the other. The modification of the state of the systemis not instantaneous, it respects the principle of causality, it is conditioned by the celerity ofthe medium. The formulation has a certain number of properties, in particular the global andlocal orthogonality of the terms in gradient and in dual curl. The discrete operators mimic someproperties of the continuous operators ∇ × ( ∇ φ ) = 0 and ∇ · ( ∇ × ψ ) = 0 , irrespective of thepolygonal geometric topology or structured or unstructured mesh.The retarded potentials are the compression and shear-rotation energies per unit of massaccumulated from a reference instant, here , up to the instant t o : φ o = − (cid:90) t o c l ∇ · V dτ ; ψ o = − (cid:90) t o c t ∇ × V dτ (7)At each instant, these retarded potentials are updated with dφ and d ψ . However, theseaccumulations are complete only in the case where the medium is perfectly elastic. Real wavespropagate with an attenuation which depends on the medium; in a Newtonian fluid the transversewaves are attenuated in a elapsed time of the order of − to − second, where they aredissipated in the form of heat. The term (cid:0) ψ o − dt c t ∇ × V (cid:1) must be replaced in this case by − ν ∇ × V , where ν is the kinematic viscosity. In the general case the discrete equation of motionthen takes the form: γ = −∇ (cid:0) φ o − dt c l ∇ · V (cid:1) + ∇ × (cid:0) ψ o − dt c t ∇ × V (cid:1) + g s (1 − α l ) φ o − c l dt ∇ · V (cid:55)−→ φ o (1 − α t ) ψ o − c t dt ∇ × V (cid:55)−→ ψ o (8)where α l and α t are the attenuation factors of the longitudinal and transverse waves. The sourceterm is written as a Helmholtz-Hodge decomposition, g s = −∇ φ s + ∇ × ψ s . The (cid:55)−→ symbolrepresents the upgrade of the potentials from time t o to time t .It is necessary to underline the autonomous character of the equation (8); indeed, the massconservation equation is not associated with the motion equation, as it is for the Navier-Stokesequation. Mass or density does not appear in equation (8), but both energy and mass areconserved. Finally, it does not contain any constitutive law.7he equivalence between mass and energy can be found through a simple analysis in orderof magnitude. If the transverse waves are very quickly attenuated in an isotropic fluid, thelongitudinal waves propagate at celerity c = c l , the speed of sound. By considering the quantity φ o − dt c l ∇ · V , it is possible to fix the order of magnitude of the second term; if V is of orderof magnitude of v , its divergence is of order v/d , and if the wave travels distance d over timeperiod dt , this term is then order of magnitude of c v and the compression energy per unit mass φ o ∝ c v . By writing the energy e = m φ o with m the mass in motion, we obtain e = m c v and,if v = c , we find the formula of the special theory of relativity e = m c for light. Solving a wavepropagation problem using the discrete equation (8) provides exactly this result.Consider the special case where c l = c t = c and apply the vector calculus formula: ∇ V = ∇∇ · V − ∇ × ∇ × V (9)By using the definition of the displacement U = U o + V dt we obtain a form which containsthe retarded potentials at the second member: d U dt − c ∇ U = −∇ φ o + ∇ × ψ o (10)The first member of this equation is a d’Alembertian (cid:3) U : c ∂ U ∂t − ∇ U (11)Thus the discrete equation of motion is, in fact, an equation of propagation of waves atvelocity c . Note, however, that c l and c t are different in the general case and the relation (9) isno longer applicable. The second Lamé coefficient ν for solids becomes the kinematic viscosityfor fluids, which does not conserve the shear energy. Mass or density is not a variable in discrete mechanics; it is replaced by the total energyper unit mass, φ o for compression energy and ψ o for shear energy. Equation (8) is autonomous:the variables ( V , φ o , ψ o ) and the celerities are sufficient to describe any problem in fluid orsolid mechanics. The conservation of mass is, however, an essential principle in mechanics; theLagrangian law of conservation, which depends only on the divergence of velocity: dρdt = − ρ ∇ · V (12)It is necessary to dissociate the variations in density due to the changes to the variables ofthe problem, here velocity, from those related to the advection of the media. This last, verydifferent, phase can be solved in multiple ways: (i) by transforming the particle derivative into atime derivative or (ii) using a posteriori one of the many Eulerian or Lagrangian methodologiesof phase transport (Volume Of Fluid, Moment Of Fluid, Level-Set, Front-Tracking, etc.).The law (12) can thus be integrated on a trajectory; the divergence of the local velocity ofthis law is also that of equation (8), updated at instant t o + dt ; explicit integration provides thesolution: ρ = ρ o e −∇ · V dt (13)Although ρ is not associated with the equation of motion, it can be computed a posteriori and be the object of a possible transport without diffusion. However, this calculation is not aresolution of the continuity equation but merely an update of the density from ∇ · V . Both8or compressible and incompressible flows, the conservation of mass is ensured intrinsically byequation (8) through that of the energy per unit of mass.Pressure p is not a retained quantity either, it is replaced by the scalar potential φ . It canalso be calculated in the form p = ρ φ in the case of ideal gases for an isothermal evolution or by p = ρ γ φ for an isentropic evolution. In general, the constitutive law p ( ρ, T ) intervenes at thislevel simply by assigning to φ a definition compatible with the behavior of the media considered.The equation of motion remains a generic kinematic law independent of any constitutive relation. Inertia is a particularly important concept in fluid mechanics because it conditions the chaoticor even turbulent evolutions of certain solutions. The physical analysis set out in reference [15]shows that it is possible to express inertia starting from a quantity named potential of Bernoulli φ B = (cid:107) V (cid:107) / in the form of a Helmhlotz-Hodge decomposition. This scalar potential is definedin all space and inertia appears in discrete mechanics as its curvature; the notion of inertia isstill complex in the current view of mechanics.In the framework of continuum mechanics, the inertial term coming from the materialderivative is written V · ∇ V or ∇ · ( V ⊗ V ) − V ∇ · V or ∇ ( (cid:107) V (cid:107) / − V × ∇ × V . Thelast term, the vector of Lamb [20], [21], is also not adapted to a discrete description of mechanics[10] since it is not a dual curl. It is therefore necessary to characterize the inertia in a differentway and to rewrite the expression of the material derivative starting from potentials φ i , theinertial scalar potential and ψ i the inertial vector potential. Physically, the variation of velocityover time on the segment Γ with respect to the material derivative, i.e. following the segmentduring its motion, is due to two actions, a direct one associated with the scalar potentials atboth extremities and another induced by the circulation of the vector potential along the contour Σ . The material derivative is fixed by the accelerations imposed from the outside and thus thevelocity variation is reduced by a quantity representing inertia. The intrinsic acceleration of thematerial medium or the particle then takes the form: γ = d V dt ≡ ∂ V ∂t + ∇ φ i − ∇ × ψ i = ∂ V ∂t + ∇ Å (cid:107) V (cid:107) ã − ∇ × Å (cid:107) V (cid:107) n ã (14)The application of the divergence and dual curl operators to the two terms of inertia makesit possible to eliminate one of these two terms this is an important difference with the mechanicsof continuous media. Equations in physics are often written as Lagrangians or Hamiltonians. Noether’s theorem[16] establishes the equivalences between the Lagrangian invariances of a system and the conservationlaws. The discrete equation can be seen as the sum of two Lagrangians which translate theexchanges between kinetic energy and potential energy like oscillators. The retarded scalarpotential φ o represents the potential energy and the term dt c l ∇ · V materializes kinetic energy;in the same way, ψ o and dt c t ∇ × V form the oscillator of the effects of rotation.The vector equation (8) is invariant to translational and rotational motions. Galilean invarianceresults in the fact that a uniform translational motion at constant velocity does not changethe equation. This property has been extended in discrete mechanics to rotational motions atconstant angular velocity [15]. The theorem of Noether serves to conclude respectively on theconservation of momentum (here the acceleration) and on that of angular momentum (here perunit of mass). The space is then qualified as homogeneous and isotropic.The invariance by translation in time expresses the fact that time passes uniformly andthat the laws of physics do not depend on it. It reflects the conservation of energy. The discrete9quation of motion is a conservation of acceleration but also of mechanical energy. The addition ofother Lagrangians on thermal or electromagnetic energy would extend the scope of this equation. The return to a discrete view of the equation of motion occurs through natural operators
DIV and
GRAD , initially described by M. Shahkov [22] for transport equations and then extendedto the curl [23]. Since then, many works of the mimetic method [24] type applied to differentequations (Navier-Stokes, Maxwell, etc.) have shown its efficiency by giving a conservativeinterpretation of existing methodologies. Mimetic methods are robust and accurate methodswhich conserve the fundamental properties of equations. At the same time, methods based ondifferential geometry and external calculus [11, 25] in a discrete form (DEC, Discrete ExteriorCalculus) have been developed over several decades. It is possible to use similar formulationsto extract the potentials and the components of a vector on polyhedral meshes by DiscreteHelmholtz-Hodge Decomposition [26, 27]. Other variants of mimetic methods have recently beendeveloped, for example for spectral [12] or finite difference [28] methods to solve the Navier-Stokesequation.The mimetic method is particularly well suited to the discrete equation of motion, written apriori as a Helmholtz-Hodge decomposition of acceleration. The numerical formulation, albeitapplying to a different physical model, is very similar to these methodologies. We will also useoperators from the mimetic methodology to transform equation (8) into a discrete equation wherethe primary operators are written G RAD , C U RL , D IV for the primal grid and ·(cid:0) G RAD , ‡ C U RL , fl D IV for the derived dual operators.The objective of the mimetic method is to create discrete approximations that preserveimportant properties of continuum equations on general polygonal and polyhedral meshes. First,we focus on the discretization of differential operators. This phase, called reduction, consists indiscretizing the continuous fields into discrete fields. The integral over a polyhedral cell is thesum of the integral over its simplexes. The reduction operator does not introduce any error inthe sense that it commutes with respect to differentiation; this result can be obtained from thegeneralized Stokes theorem. The second phase, called reconstruction, allows an approximation ofthe constitutive relations. It allows the approximate representation of dual differential operators.The design of the reconstruction operator defines the convergence rate of the numerical method.Numerical methods exist that use the mimetic framework; we can mention the mimetic finitedifference, the finite element method or the spectral element method.The primal manifold M is composed of d vertices, oriented d edges Γ , d facets of normals n and d cells. The dual mesh (cid:102) M includes d dual facets including d = 2 d in two dimensionsand d dual volumes. Using mimetic formalism, the discrete equation of motion becomes: γ = −G RAD Ä φ o − dt c l fl D IV ‹ V ä + ‡ C U RL Ä (cid:102) ψ o − dt c t C U RL V ä + g s (1 − α l ) φ o − c l dt fl D IV ‹ V (cid:55)−→ φ o (1 − α t ) (cid:102) ψ o − c t dt C U RL V (cid:55)−→ (cid:102) ψ o (15)This form causes confusion, especially on the continuous curl operator, ∇× , as we do notknow whether it applies to a vector or to a pseudo-vector. These discrete operators mimic theproperties of the continuum [24], in particular fl D IV ‡ C U RL ‹ ψ = 0 and C U RL G RAD φ = 0 . The10article derivative is expressed in the same way, with the same operators of a Helmholtz-Hodgedecomposition: γ = ∂ V ∂t + G RAD Å (cid:107) V (cid:107) ã − ‡ C U RL (cid:32) fl (cid:107) V (cid:107) (cid:33) (16)The inertial potential [15], φ i = (cid:107) V (cid:107) / , is semi-implied in time for an introduction of theinertial term within the linear system.The time-stepping procedure is close to that of classical methods: the time derivative ∂ V /∂t ≈ δ V /δt is discretized using a second-order Gear scheme. What is appreciably different, however,is the presence of time lapse dt between two mechanical balances of the system. It must bestrictly the same as the time step of the temporal discretization, δt = dt . This is one of theessential properties of the physical model presented. Discrete equation (15) makes it possibleto understand the physics of phenomena at all time scales and to capture waves at very highfrequencies, those of light for example, or to represent stationary flows. The choice of dt by theuser must simply be compatible with the physics he or she wishes to simulate. The resolution of a linear system whose unknowns are the components V of velocity is carriedout directly from the discrete system (15). The four operators of this equation, divergence,gradient, primal and dual curls, are simply formulated implicitly within the linear system. Inorder to make their description clearer, an explicit version is given. Beforehand, it is essential tonote that the physical properties, longitudinal celerity and kinematic viscosity ν which replacesthe grouping dt c t for a fluid, are constant on their respective locations, the vertex for c l and thebarycenter of the facet for ν . In this way, quantity dt c l fl D IV ‹ V is constant locally and ν C U RL V is constant on each facet.The gradient of a scalar, G RAD φ for example, is a direct operation which can be interpretedas the difference in potential at the vertices, i.e. G RAD φ = ( φ b − φ a ) /d . This quantityrepresents the component of the gradient on Γ , but not the gradient vector in space. It shouldbe remembered that the sum of accelerations is carried out on segment Γ and that at no timeis it necessary to resort to a representation in terms of a vector. Of course, the operator ∇ V ofcontinuum mechanics is not required.The primal curl at the discrete level is represented by ψ , a scalar located on facet S , obtainedfrom the circulation of the components of velocity on the primal contour Γ ∗ ; Stokes’ theoremspecifies that the curl of a vector on surface S can be calculated from its components on itscontour. Quantity ψ = ν C U RL V is thus a piecewise constant on the oriented facet. Like thediscrete gradient, the primal curl is an exact operation.The dual curl resulting from the transformation of the primal curl ψ = ν C U RL V correspondsto the reconstruction phase, which is accompanied by errors due to interpolations. This face,dual-edge transformation, leads to the dual curl ‹ ψ which allows the calculation of the circulationalong the dual δ contour. The discrete result located on the dual face is then transformed into avector carried by the primal segment Γ . It is useful to reiterate that there is no 2D/3D distinction.In two dimensions of space, for a primal planar geometry, the dual vector ‹ ψ is not in this flatplane; it is orthogonal to it, and for the planar facets in figure (1), it is directed by unit vector n . Divergence fl D IV ‹ V is calculated from velocity V associated with the segment of the primalmesh, which is transformed into dual velocity ‹ V by associating the areas of the dual facet withit. The integral of the discrete divergence on the dual cell divided by its volume makes it possibleto assign this divergence to the corresponding vertex.11n vertex-based potential-circulation, these four operators are then assembled two by two toform the two accelerations, G RAD ( fl D IV ) ‹ V ) for compression and ‡ C U RL ( C U RL V ) for rotation. G RAD ( fl D IV ‹ V ) and ‡ C U RL ( C U RL V ) In order to check the properties of each of these operators, the choice of a uniform mesh isadopted to separate the sources of intrinsic errors linked to the operators from those generatedby the deformation of the meshes.A large number of one-phase and two-phase flow test cases in fluid mechanics show thatthe convergence rate of solutions is equal to 2.0 in space and time. This is the case for theunsteady Green-Taylor vortex conducted in reference [29], but this problem is repeated innumerous publications, [30] for example. This exact time-dependent solution of the Navier-Stokes equations is used to specify the precision of the operators of the discrete formulation ona domain
Ω = [ − . , . with time steps small enough to saturate the error in time. Figure(2) shows the solution obtained for mesh of n = 16 , (i) streamlines and mesh, (ii) ψ , vectorpotential and (iii) φ scalar potential with φ = φ oB − (cid:107) V (cid:107) / where φ B is Bernoulli potential.(a) (b) (c) Figure 2.
Green-Taylor vortex on a Cartesian mesh with 256 cells, 544 edges and 289 vertices,(a) streamlines and mesh, (b) ψ = ν C U RL V , vector potential, (c) φ = φ oB − (cid:107) V (cid:107) / , scalarpotential. The scalar potential is given at the point of the primal geometric topology and the velocity isreconstructed at the center of the facets only for this graphic representation; the vector potential,also represented on the primal facets, is, on the other hand, always calculated on the barycenterof this representation and is seen like a vector, a component carried by n . The convergence studyis carried out for a spatial approximation of n = 4 to n = 1024 cells.The rate of convergence of the solution on velocity V , scalar potential φ and vector potential ψ is equal to . (figure 3). The analysis of the errors of each operator of the discrete equation(8) is carried out a posteriori from those on velocity. The scalar potential is upgraded from thedivergence of velocity φ o − dt c l fl D IV ‹ V and that of the vector potential deduced from the curl (cid:102) ψ o = − ν C U RL V . For the adopted Cartesian uniform mesh, the errors introduced by theseoperators do not modify the rate of convergence, also equal to two. Likewise, the computationof the gradient of the scalar potential and of the dual curl of ψ o saves the spatial convergenceto order two.As it is not possible to conclude on the precise error introduced by each of these operators,the velocity field of the exact solution is projected on the primal mesh, which makes it possibleto immediately extract the introduced error. It turns out that (i) the discrete gradient is exact if12 igure 3. Green-Taylor vortex, spatial convergence in the error norm L for velocity V , scalarpotential φ , vector potential ψ = ν C U RL V and gradient of potential δ = G RADφ . the scalar potential is itself exact, and that (ii) the primal curl velocity is also exact at machineprecision. Thus the second-order errors on the velocity of the numerical simulation are to beattributed to the dual curl and to the velocity divergence. Of course, the two G RAD ( fl D IV ‹ V ) and ‡ C U RL ( C U RL V ) are implicit to assemble the linear system, whose only unknown is thevelocity component V ; upgrades on φ and ψ are indeed completely explicit. Operator Propertygradient G RAD φ φ b − φ a [Γ] exact for any polygonal or polyhedral meshprimal curl C U RL V S ] (cid:88) k ∈ Γ ∗ V k d k exact for any polygonal or polyhedral meshdual curl ‡ C U RL ‹ ψ (cid:88) k ∈ δ › ψ k (cid:102) d k O ( h ) for uniform meshdivergence fl D IV ‹ V (cid:88) k ∈ ∆ › V k (cid:102) S k O ( h ) for uniform meshmimetic fl D IV ‡ C U RL ‹ ψ = 0 exactly for any polygonal or polyhedral meshmimetic C U RL G RAD φ = 0 exactly for any polygonal or polyhedral meshorthogonality local G RAD φ · ‡ C U RL ‹ ψ = 0 exactly for any polygonal or polyhedral mesh V polynomial degree ≤ G RAD ( fl D IV ‹ V ) exact for uniform mesh V polynomial degree ≤ ‡ C U RL ( C U RL V ) exact for uniform mesh Table 1.
Properties of discrete operatorsTable (1) presents the main properties of the operators of discrete mechanics. It is of coursepossible to apply the numerical methodology presented to any geometries tessellated by verydeformed meshes and, in this case, the rate of convergence of . will not be maintained. Theprecision of the method according to the quality of the mesh is not the aim of the article. Theobjective is to know the errors introduced by the processing of the jump conditions for two-phaseflows, which it is essential to separate from the errors of the method and those introduced by13ifferent considerations linked to the transport of the interfaces, calculation of the curvatures,calculation of the velocity of interfaces, etc.Finally, we show [15] that the global structure obtained has important properties [29]: (i)the vectors G RAD φ and ‡ C U RL ‹ ψ are locally orthogonal, (ii ) the discrete operators mimic theproperties of the continuous, fl D IV ‡ C U RL ‹ ψ = 0 and C U RL G RAD φ = 0 regardless of the regularfunctions φ and ‹ ψ . The numerical resolution of equations of fluid mechanics for complex physical domains or fortwo-phase flows presents difficulties of conformity of the mesh with that of the interfaces. Butthese pitfalls are also observable in problems of mechanics of one-phase fluids. The tessellationof simple surfaces or volumes, a circle for example, is not possible with elements whose faces arebased on equilateral triangles or rectangles; the irregularity of the mesh of the physical domaininevitably leads to errors, generally due to the non-planarity of the primal and dual surfaceswhere the fluxes are calculated.To obtain this property it is necessary that, for each facet of the primal and dual geometrictopologies:• the facets S of unit normals n are planar;• the segments joining the barycenter of facet S in the middle of each segment Γ areorthogonal to it, which results in the dual surfaces ∆ also being flat.where S and Γ are defined in figure (1).These conditions can only be met for regular polygonal or polyhedral meshes, such astessellations based on equilateral triangles or rectangles in two dimensions of space, or regulartetrahedra or polyhedra based on Cartesian hexahedra in three dimensions of space. This classof meshes, qualified here as uniform, is the most capable of producing excellent-quality resultsin many fluid mechanics problems.In order to overcome the problem of the conformity of the interfaces between differentenvironments or geometric domains which are sometimes difficult or even impossible to tessellatewith these uniform meshes, the proposed method consists in imposing conditions on the cells cutby an interface or a surface of the domain. There are many similar techniques designated by theterm “Immersed Boundary Methods” [31]. Other methods, corresponding to the penalization ofsome of the terms of the equation of motion, also make it possible to obtain orders of convergenceof one or two [32, 33, 34, 35, 36]. The jump conditions proposed here are specific to the discreteformulation of the equation of motion, whose physical foundations are outlined elsewhere [10, 29].The numerical methodology associated with discrete mechanics applies to non-regular polygonalor polyhedral meshes with any number of faces which do not have the property of flatness. Inthis case, the numerical solution is not free from errors, even for solutions which correspond topolynomials of degree less than or equal to two; in general the convergence of the solution is oforder two in space. The treatment of the conditions of jump in discrete mechanics differs appreciably from thatadopted for the continuous mediums. Several procedures explain these differences: (i) massor density is removed from the equation of motion except when Archimedean acceleration isactivated in an adequate potential, (ii) the effects of compression and rotation are shown by14rthogonal terms, and (iii) the jump terms conform to the derivation of all the terms of theequation, a formal Helmholtz-Hodge decomposition.The standard treatment of two-phase flows generally shows spurious currents, especially whencapillary effects are taken into account through a source term. Reviews devoted to these effects,such as that of Popinet [37], help to understand these effects due to the modeling of one or otherof the physical phenomena included. More precise analyses on the restrictions of the chosen timestep are also available in the literature, such as that of Denner [38] on surface currents that arevery sensitive to time and space steps. All these problems would require specific analysis withinthe framework of discrete mechanics. For example, if the capillary source term is not a truegradient compensated by the pressure gradient, it is not possible to obtain a stationary solutionwithout parasitic currents. Unfortunately, the defects of the curvature calculation generaterotational. b G a G * t S S xxxxxxxxxxxxxcc' b G a t S S cd' mN d (a) (b) Figure 4.
Discrete geometric topology: (a) the surface S can be composed of two domains ○ and ○ such that S = S ∪ S . The curvilinear oriented segment Σ materializes the trace on S of the interface between the two domains cutting the edges in c and c (cid:48) . (b) the interface Σ ofnormal N passes through two coplanar facets and the Γ segment at the point c ; d and d (cid:48) are thecentroids of cells and n × t the unit vector of the direct local coordinate system ( m , n , t ) . The treatment of density jumps [[ ρ ]] , viscosity [[ ν ]] or capillary pressure [[ p c ]] is not performedin the same way in discrete mechanics. The first of these is assigned to the Γ segment of theprimal geometric topology, the viscosity jump is associated with the centroids of the cells and thecapillary pressure jump is implemented by the two terms of a Hodge-Helmholtz decompositionof the capillary potential. Density jump:
The density is associated with the scalar potential and the pressure by therelation φ o = p/ρ . Like these last two quantities the density is defined on each vertex of theprimal geometry. If the density is constant and uniform in the whole domain, the pressure will beequal to the potential φ o up to a multiplicative constant. For a two-phase flow where the fluidshave constant densities ρ and ρ , only the advection of the phases will be taken into account.In the case where the density varies for other reasons, temperature or pressure variations forexample, it will have to be upgraded from the conservation of the mass: dρdt = − ρ ∇ · V (17)For the only two-phase flow cases considered here the density jump [[ ρ ]] will naturally beequal to ( ρ b − ρ a ) for all the segments cut by the interface Σ .15 iscosity jump: The formulation (8) itself sets the condition at the boundary between twofacets whose kinematic viscosities ν and ν are different and constant over them. If the velocityfields are respectively V and V , the condition at the interface is given by: ν C U RL V = ν C U RL V (18)By definition, the curl of the velocity is constant on each facet and is the same for theviscosities of the product ν C U RL V which is none other than the vector potential ‹ ψ . Thecondition (18) is actually respected locally and implicitly in all points of the domain.When a facet shown in figure (4a) is cut by the interface Σ , it is necessary to compute thevalue of the mean viscosity ν m using the properties of the theorem of Stokes by calculating thevelocity of circulation, and so of its components, along contours Γ and Γ ν m C U RL V = 1[ S ] (cid:90) Γ ν V · t dl + 1[ S ] (cid:90) Γ ν V · t dl = ν C U RL V + ν C U RL V (19)where Γ and Γ represent the contours of surfaces S and S and V and V the tangentialcomponents of the velocity on these contours. Quantity ζ m = [ S ] / [ S ] will define the partitionfactor of the phases on the facet S . The velocity flow along the interface Σ will be computedfrom the velocity on the centroid reconstructed from the components V , i.e. ( V · t Σ ) . As thecomponents of the velocity on each edge, the viscosities ν and ν on each portion of the facet S and the trace of the interface Σ inside a cell are known, it is possible to compute the averageviscosity ν m from the expression (19). It is also possible to obtain an implicit form of ν m C U RL V and integrate it into the expression (18) to express the operator ‡ C U RL ( ν C U RL V ) within thealgebraic linear system. Capillary pressure jump:
The jump condition located on edge Γ corresponding to thecapillary effects is modeled in 3D by the two source terms of a Hodge-Helmholtz decompositionof the capillary acceleration γ c : γ c = G RAD ( σ κ ξ ) − ‡ C U RL Äfl σ κ ζ ä (20)where φ c = σ κ = p c /ρ is the capillary potential with σ = γ/ρ the surface tension per unit ofmass, κ the curvature and where ξ et ζ are phase functions located respectively on the pointsand on the barycenters of the facets. When the source term is constant, for a static equilibriumfor example, it can be interpreted as well as the gradient of a scalar function or as the dual curlof a vector potential.Like the scalar potential of the acceleration, the capillary potential is expressed in m s − .The capillary potential difference between the two fluids ∆ φ c is independent of the density, butthis is not the case for the capillary pressure which must be computed by ∆ p c = ρ ∆ φ c . Thecapillary pressure jump p c is defined from quantities located at the points of the primal geometrictopology, the potential φ c and the density ρ . This model belongs to the class of sharp interfaces.To carry out a simulation of a two-phase flow it is advisable to redefine the scalar potential φ o = p/ρ at each iteration in time in order to take account of the advection of the phases.Consider the 2D-case of a circle of curvature κ = 1 /R placed in the center of a square withside R ; the constant surface tension is equal to σ . When the discrete equation of motion isrestricted to only capillary effects it is written: γ = −G RAD Ä φ o − dt c l fl D IV ‹ V ä + G RAD ( σ κ ξ ) (21)This equation is solved from zero values of the velocity field and the scalar potential. Thesolution is obtained in a single step with a precision of ≈ − compared to the theoreticalsolution φ oc = σ κ . At the end of the resolution the velocity is zero and of course fl D IV ‹ V = 0 .16 mplementation: The general numerical methodology associated with solving two-phaseproblems whose interfaces Σ do not conform to the structured or unstructured meshes used isdescribed below:• locate the points of the primal mesh that are inside the object (phase ) by a ray-tracingmethod;• detect all intersections of the trace of the interface with Γ segments and find the intersectionpoints for each facet exactly S ;• define the polygons with 3, 4, 5 sides (or more depending on the tessellation) and calculatetheir areas; check whether the total area of the object is found exactly;• compute the occupancy rate of the ○ and ○ phases on each facet;• compute the circulation on the contour of each cell portion crossing the interface Σ ;• compute the values of the properties on each intersected element, the density ρ on thevertices, the viscosity ν m on the facet and the mean curvature on the segment.All these operations are performed with very efficient differential geometry routines. Theywill be repeated at each time step for a real simulation where the interface is moving within a fixedmesh. In the case of a Lagrangian approach (ALE for example), the interface is in conformity withthe mesh and the jump conditions are not necessary, but the formulation applies nevertheless. Standard test cases relating to one-phase or two-phase flows have already been carried out[34], with a similar formulation to that presented, in particular, on the benchmark of a bubblerising in a liquid [39]. Others have been carried out more recently for two-phase flows or fluid-structure interactions [10, 29]. The results are in agreement with those of the literature.The numerical solutions shown here generally have a solution of degree less than or equal totwo on the velocity which corresponds to φ and ψ potentials linear or piecewise linear and thusto accelerations G RAD φ and ‡ C U RL ‹ ψ constant. When each dual ∆ facet is orthogonal to the Γ segment, the numerical solution must be exact. Consider two fluids of density ρ and ρ and kinematic viscosity ν and ν ; the two fluidsare superimposed in a cavity of height H = h + h where h and h are the heights occupiedby each of the fluids. At the initial moment, the set is subjected to gravity according to y , g = − g k = − g e y where vector k is vertical ascending. The heavier fluid can logically be below,but even if the equilibrium is potentially unstable, the theoretical and numerical solutions areperfectly stable. Initially, the pressure and velocity are zero. Figure (5) shows the scalar potential φ o ( y ) solved and the pressure p ( y ) which is extracted after the explicit upgrade defined by thecontinuous expression (22): p b = p a − (cid:90) ca ρ ∇ φ · t dl − (cid:90) bc ρ ∇ φ · t dl (22)The system equation (8) is resolved with a time step of dt = 10 s in order to obtain thesteady solution in a single step. During the resolution of the linear system, the velocity is directed17 igure 5. Equilibrium of two fluids subjected to the influence of gravity; on the left, the evolutionof the potential φ o ( y ) and on the right the evolution of the pressure p ( y ) . The theoretical solutionis represented by lines and the numerical solution by points; the solution is accurate to the machineprecision. downwards for the two fluids, but as the incompressibility is applied very strongly, it is null atthe end of the algebraic resolution and the solution of problem ( φ o , V ) at equilibrium is obtainedinstantly with V = 0 ; the uniform mesh is composed of cells. The heights of the two fluids arerespectively h = 0 . m and h = 0 . m ; the heavy fluid placed in the lower part of the cavityhas a density equal to ρ = 10 and the light fluid a density of ρ = 1 .The velocity obtained is at zero divergence ( fl D IV ‹ V ≈ − ) and the pressure correspondsto the theoretical solution sought. Viscosities ν and ν are taken into account only duringresolution, as the stationary solution is independent of this physical parameter. The densityratio is chosen as being equal to ρ /ρ = 10 and may be arbitrary.This test case has already been carried out recently using a finite element method [40], wherethe analytical solution is found exactly. One of the advantages of the incompressible formulationpresented here lies in the absence of density in the equation of motion, which makes it possibleto obtain a solution on the potential φ = p/ρ regardless of this. The same test case, whose meshdoes not conform with the interface, is carried out in a square [ − . − . inclined with aslope of tg ( θ ) = 0 . , that is, an angle θ = 0 . rd and a descending vertical gravity g = − e y .The cavity is filled half with a fluid of density ρ = 10 and half of density ρ = 1 , the interfaceis horizontal and the pressure is zero initially. The flow is supposed to be incompressible andthe divergence is zero throughout the calculation. The solution can be obtained in an iterationin time or a time step dt ; the stationary solution is the same, it is determined with machineprecision whatever the adopted uniform mesh. Figure (6) shows the numerical solution obtainedfor a Cartesian mesh of cells.This case does not present parasitic currents at any time; the interface remains flat duringthe simulation. When the surface is initially inclined with respect to the horizontal direction, theArchimedes effects generate a movement which attenuates if the viscosity is present to give theprevious solution again. This dynamic case of sloshing will not be treated within the frameworkof this analysis devoted to the treatment of the conditions of jumps, but the solution obtainedwith the same model can be found elsewhere [10], where the oscillation frequencies are comparedto the results of linear theory. It should be noted that the density is not required for this problem,because the equation of the statics of fluids in a continuous medium is written −∇ p + ρ g = 0 ,while the discrete equation is written −G RAD φ + g = 0 . The quantity φ is the potential ofacceleration, which is not the case of pressure. 18 igure 6. Equilibrium of two fluids subjected to the influence of gravity; on the left the scalarpotential φ ∈ [ ± . and on the right the pressure p ∈ [0 , ; the numerical solution isequal to the theoretical solution to machine precision and the velocity is zero; the isovalues arenot drawn in the vicinity of the interface because the graphic interpolation between two points islinear and not, linear by pieces. In fact, the analytic treatment of this problem does not require any resolution, it is sufficientto note that the vector g is constant and that it can derive from a vector potential ψ g or froma scalar potential φ g such as g = G RAD φ g . We can explicitly calculate the potential on eachpoint of the structured or unstructured mesh: φ b = φ a − (cid:90) ba g · t dl (23)It is sufficient to choose a first point whose value will be arbitrary and to traverse the wholeof the primal mesh point after point to calculate the potential to a constant and the pressure p = ρ φ . The flow of Poiseuille in a planar channel can be simulated using a Cartesian mesh conformingto the geometry of the channel. In this case, the numerical solution is accurate to machineprecision. The test case presented here corresponds to the same planar channel but is inclinedin some way in a regular Cartesian mesh. The height channel h = [ ± . is bounded by twosolid planes inclined at θ = 26 . from the horizontal. The fluid has a density equal to ρ and akinematic viscosity equal to ν . The flow is stationary and incompressible.The jump condition (19) is applied in this case for cells cut by the solid wall representedby viscosity ν → ∞ defining the channel, in practice ν = 10 . A simple calculation leadsto a viscosity value ν m which depends only on the viscosity of the fluid moving in the channel, ν m = S / S where S is the portion of surface occupied by this fluid. The equation of motion isderived from discrete mechanics, but since the problem is stationary, the only physical parametersare r = dt c l and the kinematic viscosity ν = dt c t . As inertia is zero for this case, the solutioncan be obtained in a single iteration. γ = −G RAD Ä φ o − r fl D IV ‹ V ä + ‡ C U RL Ä (cid:102) ψ o − ν C U RL V ä + g (cid:102) ψ o − ν C U RL V (cid:55)−→ (cid:102) ψ o (24)19ith g = 1 e x = G RADφ g = G RAD ( − x ) = − ‡ C U RL (cid:102) ψ g = − ‡ C U RL ( y ) . The numerical parameter r = dt c l makes it possible to maintain a divergence that is as small as necessary, an order ofmagnitude fl D IV ‹ V ≈ /r . Figure 7.
Poiseuille flow in a planar channel h = [ ± . inclined to θ = 26 . ° in adomain of dimension ( ± ; potential vector of acceleration (cid:102) ψ o ∈ [ ± . ; the mesh is regularCartesian cells. The scalar potential φ o and the divergence of the velocity fl D IV ‹ V are null tomachine precision. n Residual V Residual ψ o .
30 10 − .
55 10 − .
55 10 − .
04 10 − .
51 10 − .
82 10 − Table 2.
Poiseuille flow in an inclined channel; residual on the velocity in norm L betweenthe theoretical solution projected on the segments of the primal mesh and the numerical solution;residue on the vector potential (cid:102) ψ o . The solution shown in figure (7) corresponds to an equal steady-state potential at (cid:102) ψ o = ν C U RL V and a velocity equal to V · e u = y + a y + b where e u is the axial vector of theinclined channel. The figure schematically represents the parabolic profile of the velocity and ofthe current lines between the two solid walls.Table (2) gives the residuals between the theoretical solution and the numerical simulationson the velocity and the vector potential for three different meshes, thus showing the independenceof the spatial approximation for the problem posed. For this example, the two fluids have finite viscosity of ν and ν . Since the variations inpotentials are linear, the homogenized viscosity has a simple form: ν m = Å ν S S + 1 ν S S ã − (25)As the geometry of the channel is fixed in time, it is easy to calculate ν m for each cell. In steadystate, the density does not affect the solution. 20onsider a duct of height h = 1 bounded by two horizontal surfaces of zero velocity. Thelower kinematic viscosity fluid ν = 1 and the higher viscosity fluid ν = 10 occupy the spacerespectively between y ∈ [0 , y ] and y ∈ [ y , . The fluid is set in motion by a pressure gradientor vector potential gradient in discrete mechanics, such that ‡ C U RL (cid:102) ψ o = 1 . . The interface ispositioned at value y = 0 . and does not correspond to a particular value associated with themesh. The initial velocity and scalar potential values are zero. The solution can be obtained ina single iteration since inertia is zero for the problem posed. Figure 8.
Two-phase Poiseuille flow of height h = 1 where the interface is positioned in y = 0 . .The viscosity of the fluid at the bottom is ν = 1 and that of the higher fluid is ν = 10 ; the meshis cells. Comparison between the analytical solution (line) and the numerical values (points);the divergence of the velocity fl D IV ‹ V and the scalar potential φ o are null to machine precision. The condition at the boundary between the two fluids, whose velocity fields are V and V ,is given implicitly by the discrete formulation of the equation of motion. The residue in norm L on the velocity and on the vector potential is given by table (3): n Residual V Residual (cid:102) ψ o .
45 10 − .
27 10 − Table 3.
Two-phase Poiseuille flow in planar channel; residual on the velocity in norm L between the theoretical solution projected to the barycenters of the primal mesh and the numericalsolution and residue on the vector potential (cid:102) ψ o . The solution is accurate to machine precision for all quantities, regardless of the number ofcells adopted. The test case of the horizontal planar channel is taken up in a mesh based ontriangles. With triangular meshes made with any software, several imperfections usually appear;for example the orthogonality of the segments of the primal topology and those connecting thebarycenters is not assured. If the objective is to achieve machine precision on the numericalsolution, then it is necessary that all the metrics be evaluated with the machine precision, forexample the position of the points must be given with significant digits in double precision.Generally, the results are in the order of two, but the absolute error in the vicinity of theinterface between the two fluids is very important, whatever the interpolation adopted for the21iscosity. Since the theoretical solution corresponds to two polynomials of degree equal to two,it is possible to find an exact numerical solution (figure 9).The height channel h = 1 is filled with two fluids of viscosity ν and ν with an interface Σ located at y = 0 . . The equilateral triangle mesh contains × cells. Figure 9.
Two-phase Poiseuille flow of height h = 1 where the interface is positioned in y = 0 . .The viscosity of the fluid at the bottom is ν = 1 and that of the higher fluid is ν = 10 ; the meshis × cells. Vector potential field (cid:102) ψ o and some current lines; the values of the divergence ofthe velocity fl D IV ‹ V and the scalar potential φ o are null to the machine precision. The numerical solution is represented in figure (9) which shows the mesh used, the vectorpotential field and some current lines. The interface Σ is set to y = 0 . and the stream lines areunaffected by changes in viscosity. n Residual V Residual (cid:102) ψ o .
02 10 − .
77 10 − Table 4.
Two-phase Poiseuille flow in planar channel for a uniform mesh based triangles;residual on the velocity in norm L between the theoretical solution projected to the barycentersof the primal mesh and the numerical solution and residuals on the vector potential (cid:102) ψ o . Table (4) gives the residue of the numerical solution on the vector potential and velocity. Theviscosity ratio here is equal to but it is possible to set ratios of (for a non-conforming freesurface V · n = 0 ) to infinity (for a non-conforming solid wall V = 0 ). This formulation makesit possible to process 3D interfaces without modification to address complex geometry problemswith non-conforming meshes. The Hadamard-Rybczynski solution [41, 36] corresponds to an incompressible, stationarynon-inertial flow of a bubble in a fluid of different density under the effect of the Archimedeanforce. In fact, as the densities are constant in each of the fluids, the source term to be integratedinto the equation of motion is constant. Since the external solution is in /r , the numericalsolution can only be at the order of the scheme and the absolute error then depends on it. Onthe other hand, the internal solution corresponds to a polynomial of degree equal to two andit is then possible to obtain the theoretical solution to machine precision whatever the spatialapproximation. 22he solution is sought from the equation of motion (8) solved in the reference frame of thesphere. The separation of the variables was used by Hadamard and Rybczynski to obtain thesolution: V = K (cid:2)(cid:0) r (cid:1) cos θ e r − (cid:0) r (cid:1) sin θ e θ (cid:3) ψ = − K r sin θ e ϕ (26)with r = √ x + z ( e r = e x + e z ) and K = g R ( ρ e − ρ ) / µ e +3 µ/ where R is the radius of thesphere, µ and ρ the viscosity and density of the inner fluid and µ e and ρ e those of the exterior fluid.For the simulation, the density and viscosity dependent parameter is chosen equal to K = 1 . Thesource term of the system (8) is written h = −G RAD φ g + ‡ C U RL (cid:102) ψ g = −∇ ( − r ) + ∇ × ( y e ϕ ) . The methodology used is based on an unstructured mesh notconforming to the unit circle representing the spherical bubble (treated axisymmetrically). Theequilateral triangle mesh requires the triangles cut by the circular interface Σ to be treated andthe non-intersected outer triangles to be excluded. A series of differential geometry algorithmsallows calculation of the area of polygon portions with one side curved. The reconstitution ofthe unit circle makes it possible to calculate its area A = π R to machine precision. From theinformation on the proportion of the area of the triangle belonging to the domain, it is possibleto calculate the viscosity of each of the cells cut by the interface. Figure 10.
Hadamard-Rybczynski’s solution inside the sphere; left vertical velocity for a mesh of n = 96 elements; center and right, vector potential ψ and current lines for a mesh of n = 22608 equilateral triangles. The solution is exact and the divergence is zero. The source term and incompressibility are sufficient to define this test case completely. Theseparation of the variables used by Hadamard to obtain the solution is proof that the effects ofcompression and shear are disjoint. Figure (10) shows the vertical velocity in the bubble and thevector potential ψ , exact to machine precision regardless of the spatial approximation. Valuesoutside the circle have no meaning. Since the current function is not a primitive variable ofdiscrete mechanics, it is recalculated on the points by solving a Poisson equation; the currentfunction is indeed only one component of the vector potential of the velocity, whereas ψ is thevector potential of the acceleration. The velocity components and vector potential ψ obtained bynumerical simulation are accurate to machine precision regardless of the spatial approximation.The solution is represented in figure (10). It is obtained in a single step and a few conjugategradient iterations. 23 D-solution:
The solution to this problem in a three-dimensional space is treated using astructured mesh based on regular hexahedrons, where the outer facets not concerned by the flowin the bubble are eliminated it a priori. As regards the principle of the methodology, nothingis changed from the 2D-axisymmetric approach and the same procedure is applied whatever theorientation of the facets. The trace of the sphere of unit radius on the facets makes it possibleto calculate the areas S of the parts of the spherical domain. The boundary conditions on thesphere naturally correspond to a slip of the fluid. Figure 11.
Hadamard-Rybczynski 3D solution inside the sphere; the mesh is composed of hexahedral cells. Stream lines and equipotentials show that the velocity field and the vectorpotential ψ are axi-symmetric. The stationary solution is exact and the divergence is zero. Figure (11) represents the solution obtained directly with system (8) graphically. The solutionobtained on the velocity field is exact and the vector potential ψ is strictly linear per surface.This example treated with a mesh conforming to the unstructured sphere inevitably leads toerrors due to the tessellation of the sphere surface and the polyhedral volume mesh. The choiceof a Cartesian structured mesh is of considerable interest for various reasons and the applicationof these jump treatment techniques makes it possible to solve one of its difficulties, consisting intaking account of the complex geometry. The incompressible equilibrium of a droplet subjected to capillary acceleration has absolutelyno dependence on the viscosity or density of the fluids. The pressure difference between fluidsis given by Laplace’s law ∆ p = γ κ , or, in discrete mechanics, by ∆ φ o = σ κ . Indeed, anincompressible movement, linked to an infinite celerity of sound c → ∞ , leads to an instantaneousequilibrium. The equation (8) in the presence of a capillary force differs significantly from theNavier-Stokes equation; in particular the conventional couplings between velocity and pressureare replaced by an integration of the incompressibility constraint in the equation itself. Thisproblem does not generate a priori any rotational one usually interpreted as spurious currents.The equation of motion is in this case: γ = −G RAD Ä φ o − r fl D IV ‹ V − σ κ ξ ä (27)The steady solution of this problem is obtained explicitly by considering the medium asincompressible reducing the equation (27) to −G RAD ( φ o − σ κ ξ ) = 0 ; the two terms being24radients the scalar potential of the acceleration can be obtained explicitly by calculating thesolution, up to a constant, from one point to another by the relation: φ ob = φ oa − (cid:90) ba ∇ ( σ κ ξ ) · t dl (28)where ξ is the phase function equal zero or one. This solution is represented in figure (12)for three meshes, a Cartesian mesh, another an unstructured mesh of triangles and a third onequilateral triangles. Figure 12.
Equilibrium of a droplet under a capillary effect simulated using meshes based onCartesian quadrangles, triangles and equilateral triangles; the potential difference on either sideof the interface is exact.
The simulation with all terms of equation (8) is performed from null initial pressure andvelocity fields. When the quantity dt c is very large, the divergence becomes very small and thesolution of −G RADφ o + G RAD ( σ κ ξ ) = 0 is obtained instantly. A time step of dt = 10 s allowsthe steady solution to be reached instantaneously where V = 0 and φ o = σκ . The accuracy ofthe solution is order of magnitude of machine precision; the difference between analytical andnumerical solution it is kept smaller than − after 2 iterations. The physical model presentedis suitable for solving various problems in mechanics. For a fixed intrinsic celerity, the natureof the flow depends on the observation time scale dt ; for example, the movements of water, anessentially incompressible liquid ( c = 1500 m s − ) , cause acoustic waves to appear with weaktime constants. Excluding microscopic interactions, the model (8) is representative of physicalphenomena at all scales of time and space.The model (8) is full-compressible and the behavior of the solutions is closely associatedwith the longitudinal and transverse waves which propagate in the media. For the two fluidsof the Laplace test case, the transverse waves attenuate very quickly but there remain thelongitudinal waves which propagate at the celerity of sound c = 1 / √ ρ χ T , where χ T is thecompressibility coefficient. Thus, if the scalar potential is zero at the initial moment, the capillaryacceleration G RAD ( σ κ ξ ) causes a radial collapse of the drop on itself, compensated by thequasi-incompressibility of the fluid; as the fluids have finite celerities the stationary solution isnot reached instantly if the time steps dt are lower than R/c . Figure (13) shows the solution attime t = 5 10 − s obtained with a time-step dt = 10 − s for two non-viscous fluids.As the two fluids have different properties, the radial expansion wave in the external fluidand the compression wave in the drop propagate at different celerities, since the square cavityis closed. Continuing the simulation leads to a succession of compressions and detents whichinteract with the walls of the cavity. For inviscid fluids, the kinetic energy should be conservedover time provided that the phenomenon is simulated with sufficient temporal precision.25 igure 13. Droplet of radius R = 2 10 − m in a square cavity filled with two inviscid fluids;scalar potential field; Cartesian mesh of cells, dt = 10 − s , t = 5 10 − s with c = 10 m s − outside and c = √ m s − inside; the outline of the drop is represented by a white circle. Figure 14.
Kinetic energy in a droplet of radius R = 2 10 − m in a square cavity filled withinviscid fluids; norm of velocity, mesh of and time-step dt = 10 − s , with various celeritiesof sound (a) c = 10 m s − , (b) c = 10 m s − , (c) c = 10 m s − Depending on the values of the celerity of sound, the steady solution of capillary equilibriumis obtained more or less quickly. The oscillations of the velocity norm represented in figure (14)are absolutely not parasitic currents but reflect the temporal variations of the potential φ ofthe compressible movement. The incompressible motion c → ∞ must lead to an instantaneouscapillary equilibrium and the incompressible model causes the instantaneous extinction of allacoustic waves.We can see that (i) the solution does not depend on the main mesh, (ii) the curvature26btained with differential geometry routine is exact , (iii) the oscillations of the velocity normare independent of the viscosity, (iv) the solution is independent of density.The interpretation of the results of the dynamic and static equilibrium of a drop subjectedto a capillary pressure jump differs from that generally given for simulations carried out in thecontext of a continuous medium for an incompressible movement. The Capillary Ca , Laplace La or Ohnesorge Oh numbers representative of this problem involve the viscosity of fluids and theirdensity.The Continuum Surface Force (CSF) model generally introduced in continuum mechanicsto describe the capillary term as γ κ δ Σ n = γ κ ∇ c has no reason to be a gradient, in general.The spurious currents observed by many authors are due to curvature calculation errors whichgenerate vorticity whose intensity depends on the numerical methodology used. The oscillations of a drop subjected to capillary accelerations are an opportunity to highlightthe originality which consists in removing the density of the quantities present in the incompressibleNavier-Stokes equation. Indeed, the two physical quantities present in the discrete equation (8)are the celerities c l = 1 / √ ρ χ T and c t = (cid:112) µ/ρ which show three physical quantities, ρ , µ as wellas χ T whose measurement is generally carried out indirectly by that of the celerity. For timeconstants greater than τ ≈ − s , the grouping dtc t must be replaced by the kinematic viscosity ν for homogeneous and isotropic fluids. In addition to the reduction of physical quantities, onecan observe with table (5) that, in the case of a two-phase water-air flow, the ratios of the selectedquantities are much lower than those adopted conventionally. ρ χ T µ c l √ ν water − − − air − −
316 3 .
16 10 − ratio − − .
16 0 . Table 5.
Physical parameters of two fluids close to those of water and air.
The historical work of Lord Rayleigh [42, 43] on the oscillations of drops and bubbles served todefine the frequency of infinitesimal-amplitude oscillations of incompressible liquid drops aboutthe spherical shape in vacuum. More recent studies, including [20], [44], [45], [46], have extendedthe field of investigation to dissipative effects using linear theories. The problem of the oscillationof a droplet in nonlinear theory for a small Ohnesorge number has been addressed recently[47] using the potential flow assumption to reduce the corresponding free boundary problemformulated on a time-dependent domain into a nonlinear system of integro-differential equations.For two-dimensional planar drops oscillating about a circle, the frequency of the oscillations isgiven by: f th = m ( m − γρ R (29)where m is the mode number and the density ρ is that of the drop fluid. The frequency isindependent of the viscosity but the attenuation of the oscillations depends on it.The test case is chosen so as to show that the formulation of the jumps in density, viscosity andcapillary pressure, as well as the numerical methodology implemented, do not affect the propertiesof spatial convergence of the system (8). Consider a circular drop of radius R = 210 − m deformedinto an ellipse with a radius ratio equal to R x /R y = 1 . at the center of a square of dimension27 = 10 − m . The properties of the drop recalled in table (5) are very close to those of water,and the physical characteristics of the gas which surrounds it correspond substantially to thoseof air. The surface tension fixed at γ = 2 10 − kg s − or σ = 2 10 − m s − is chosen small enoughto avoid confusing the capillary effects and the acoustic waves of the oscillations of the ellipseover time.The simulations are performed on Cartesian meshes with degrees of freedom varying from to ; the initial ellipse is represented by a string of , or markers. Figure (15)shows the fields of scalar potential φ o and instantaneous vector potential ψ o for a time close tothe initial state. Figure 15.
Snapshot of scalar potential φ o , vector potential ψ o and some trajectories for a time t = 10 − s and a spatial approximation n = 512 The variations of the curvature within the source term G RAD ( σ κ ξ ) generate an accelerationand a movement which tends to reduce the curvature, but the inertia leads to oscillations whichslowly dampen given the low viscosity of water. This case test makes it possible to validate,among others, the specific nonlinear terms of the discrete formulation. An example of theevolution of kinetic energy in the field is shown in figure (16). Figure 16.
Evolution of the kinetic energy E k over time of an elliptical drop in a squarecontaining gas for σ = 2 10 − m s − , ρ = 10 kg m − , dt = 10 − s and spatial approximation n = 64 .
28e can clearly observe the two behaviors of the movement of the drop over time, first inertialand then viscous; the two lines in E k ∝ exp ( − α t ) represent their respective attenuations. Formuch greater times the velocity tends towards zero without parasitic currents and the difference inscalar potential takes the theoretical value ∆ φ o = 0 . to exactly satisfy the equality −G RADφ o + G RAD ( σ κ ) = 0 .The simulations carried out, with sufficiently weak time steps (10 − s ) to neglect the errors intime and different spatial approximations, make it possible to calculate the rate of convergenceof the numerical solution; this is given by table (6) and figure (17). n f ∞ orderfrequency .
928 36 .
880 39 .
154 40 .
132 40 .
560 40 .
669 40 .
706 1 . Table 6.
Frequency of oscillations f with spatial approximation n . The rate of convergence in space calculated from the value f ∞ , obtained using Richardson’sextrapolation, is close to ; the theoretical value (29) for the mode for m = 2 gives the value of f th = 38 . s − . This difference can be explained by the nature of the problem treated, a littledifferent from the non-viscous theory, the marked nonlinear effects, the confinement due to thepresence of walls, and of course the numerical errors of discretization of the interface. Figure 17.
Convergence of oscillation frequency f with spatial approximation n . This example is not a complete study of the free oscillations of drops, but it does serve toverify that the numerical methodology presented simulates the phenomenon with an acceptableprecision. However, it is possible to affirm that the initial objective has been reached: multipleerrors due to the processing of jump conditions, to the transport of interfaces by a front-trackingmethod, to the calculation of curvatures, to the compressible formulation, etc., are lower thanthe precision of the discrete formulation itself.
In order to give an idea of the robustness of the formulation, an example of two-phase flow isperformed from the real values of the physical properties, density, viscosity and surface tension29hich vary significantly, spatially and temporally.(a) (b)
Figure 18. (a) Shape of the bubble at different times: a : 0 s , b : 0 . s , c : 1 s , d : 1 . s , e : 2 s , f : 2 . s , g : 3 s , h : 3 . s , i : 8 s ; (b) snapshot of the vertical component of velocity at t = 0 . s ,trajectories and shape of the bubble extracted from the position of the markers. A closed cavity of height H = 0 . m and width L = 0 . m is filled with mercury; an airbubble of radius R = 0 . m is initially positioned in ( x = 0 . , y = 0 . within the mercuryand the whole is subjected to vertical gravity of g = − · e y ms − . The properties of mercury (inSI units) are as follows: dynamic viscosity µ = 0 . , density ρ = 13600 , longitudinal celerity c = 1407 . Those of air are dynamic viscosity µ = 1 .
85 10 − , density (constant) ρ = 1 . ,longitudinal celerity c = 293 . The surface tension between the two fluids is taken as γ = 0 . .The boundary conditions of the cavity are adherent V = 0 . The values are given for informationonly and may vary according to the chosen methodologies. The time step is dt = 10 − s andthe total time of simulation is s . The regular Cartesian mesh adopted is · cells. Theair bubble is initialized by a string of markers defining the radius circle R . The position ofthe interface on Γ segments is given by fast routines of differential geometry; it allows precisecalculation of the density on each segment. The curvature is calculated from the dihedral anglebetween three consecutive markers.The evolution of the shape of the bubble in time is shown by figure (18a) at different times.Figure (18b) shows a snapshot of the shape of the bubble at t = 0 . s , some trajectories andthe vertical component of the velocity field. The physical parameters lead to a divergence ofless than − throughout the calculation. The flow behind the bubble is highly inertial andunsteady. Even though the bubble motions weaken strongly after s in the vicinity of the upperwall, the vortex flows persist for a very long time, given the kinematic viscosity of mercury ( ν ≈ − m s − ) . The field of the scalar potential is of order one and perfectly continuous; itallows a return to the pressure, only if necessary.30 Conclusions
The first part of this article presents the physical model, called discrete mechanics, andthe numerical methodology which is closely associated with it. The general principles of thederivation of the law of motion abandon the notion of continuous medium to directly constructthe discrete formal framework to establish a vector equation expressing the fact that the intrinsicacceleration of a material medium is equal to the accelerations applied to it in a given direction.This alternative formulation to the Navier-Stokes equations makes it possible to find the resultsof these with constant physical properties. The application of this new formulation to two-phaseflows requires the establishment of conditions of jumps of the physical properties in line with thechoice of expressing each physical effect in two contributions, one with curl-free and the secondwith divergence-free.The formulation is mainly applied to very simple test cases where a theoretical solution ofdegree two is available. The resolution of these two-phase cases makes it possible to restore thetheoretical solution with errors of the order of magnitude of the machine accuracy (without anynumerical artifact), regardless of the chosen spatial approximation. The treatment of contactdiscontinuities does not affect the accuracy of the basic scheme. In the other cases [10, 29] andin general, this formulation is of order two in space and time. Its robustness can be evaluatedon a flow with a high contrast of physical properties.
Author Contribution
Author: Physical modeling, Conceptualization, Methodology, Research code, Validation,Writing- Original draft preparation, Reviewing and Editing.The paper has been checked by a proofreader of English origin.
Declaration of competing interest
There are no conflict of interest in this work.
References [1] R. Scardovelli, S. Zaleski, Direct numerical simulation of free-surface and interfacial flow,Annual Review Fluid Mech. 31 (1999) 567. doi:10.1146/annurev.fluid.31.1.567 .[2] A. Prosperetti, G. Tryggvason, Computational Methods for Multiphase Flows, CambridgeUniversity Press, Cambridge UK, 2007.[3] R. Fedkiw, T. Aslam, B. Merriman, S. Osher, A non-oscillatory Eulerian approach tointerfaces in multimaterial flows (the ghost fluid method), J. Comput. Phys. 152 (1999)457–492.[4] C. Wang, T. Liu, B. Khoo, A real ghost fluid method for the simulation of multimediumcompressible flow, SIAM J. Sci. Comput. 28 (2006) 278–302.[5] A. Guittet, M. Lepilliez, S. Tanguy, F. Gibou, Solving elliptic problems with discontinuitieson irregular domains - the Voronoi interface method, J. Comput. Phys. 298 (2015) 747–765.[6] C. Liu, C. Hu, A second order ghost fluid method for an interface problem of the poissonequation, Commun. Comput. Phys. 22 (2017) 965–996.[7] P. Trontin, S. Vincent, J. Estivalezes, J. Caltagirone, A subgrid computation of the curvatureby a particle/level-set method. Application to a front-tracking/ghost-fluid method for31ncompressible flows, J. Comput. Phys. 231 (20) (2012) 6990–7010. doi:10.1016/j.jcp.2012.07.002 .[8] P. P.G. Tucker, Z. Pan, A cartesian cut cell method for incompressible viscous flow, AppliedMathematical Modelling 24 (2000) 591–606. doi:10.1016/S0307-904X(00)00005-6 .[9] D. Ingram, D. Causon, C. Mingham, Developments in Cartesian cut cell methods,Mathematics and Computers in Simulation 61 (2003) 561–572. doi:10.1016/S0378-4754(02)00107-6 .[10] J.-P. Caltagirone, Discrete Mechanics, concepts and applications, ISTE, John Wiley & Sons,London, 2019. doi:10.1002/9781119482826 .[11] M. Desbrun, A. Hirani, M. Leok, J. Marsden, Discrete exterior calculus,arXiv/math/0508341v2 (2005) 1–53.[12] A. Palha, M. Gerritsma, A mass, energy, enstrophy and vorticity conserving (MEEVC)mimetic spectral element discretization for the 2D incompressible Navier-Stokes equations,J. Comput. Physics 328 (2017) 200–220. doi:10.1016/j.jcp.2016.10.009 .[13] J. Maxwell, A dynamical theory of the electromagnetic field, Philosophical Transactions ofthe Royal Society of London 155 (1865) 459–512.[14] C. Will, Theory and Experiment in Gravitational Physics, Cambridge University Press,Cambridge, United Kingdom, 2018.[15] J.-P. Caltagirone, On Helmholtz-Hodge decomposition of inertia on a discrete local frameof reference, Phys. Fluids 32 (2020) 083604. doi:10.1063/5.0015837 .[16] Y. Kosmann-Schwarzbach, Nother Theorems. Invariance and Conservations Lawsin the Twentieth Century, Springer-Verlag, New York, 2011. doi:10.1007/978-0-387-87868-3 .[17] J. Guermond, P. Minev, J. Shen, An overview of projection methods for incompressibleflows, Comput. Methods Appl. Mech. Engrg. 195 (2006) 6011–6045. doi:10.1016/j.cma.2005.10.010 .[18] A. Liénard, Champ électrique et magnétique produit par une charge électrique concentréeen un point et animée d’un mouvement quelconque, L’Éclairage électrique 27 (1898) 5–112.[19] J.-P. Caltagirone, Discrete Mechanics, ISTE, John Wiley & Sons, London, 2015. doi:10.1002/9781119058588 .[20] H. Lamb, Hydrodynamics, e edition, Dover, New-York, 1994.[21] C. Hamman, J. Klewick, R. Kirby, On the Lamb vector divergence in Navier-Stokes flows,J. Fluid Mech. 610 (2008) 261–284. doi:10.1017/S0022112008002760 .[22] M. Shaskov, Conservative Finite-Difference Methods on General Grids, Boca Raton: CRCPress, 1996. doi:10.1201/9781315140209 .[23] J. Hyman, M. Shashkov, Natural discretizations for the divergence, gradient ans curl onlogically rectangular grids, SIAM J. Num. Anal. 36 (1999) 788–818.[24] K. Lipnikov, G. Manzini, M. Shashkov, Mimetic finite difference method, J. Comp. Physics257 (2014) 1163–1227. doi:10.1016/j.jcp.2013.07.031 .3225] M. Meyer, M. Desbrun, P. Schröder, A. Barr, Discrete Differential-Geometry Operators forTriangulated 2-Manifolds, In: Hege HC., Polthier K. (eds) Visualization and MathematicsIII. Mathematics and Visualization, Springer, Berlin, Heidelberg, 2003.[26] E. Ahusborde, M. Azaiez, J.-P. Caltagirone, A primal formulation for the Helmholtzdecomposition, J. Comp. Physics 225 (1) (2007) 13 – 19. doi:10.1016/j.jcp.2007.04.002 .[27] A. Lemoine, J.-P. Caltagirone, M. Azaiez, S. Vincent, Discrete helmholtz-hodgedecomposition on polyhedral meshes using compatible discrete operators, Journal ofScientific Computing 65 (2015) 34–53. doi:10.1007/s10915-014-9952-8 .[28] A. Abbà, L. Bonaventura, A mimetic finite difference discretization for the incompressibleNavier-Stokes equations, International Journal for Numerical Methods in Fluids 56 (2008)1101–1106. doi:10.1002/fld.1678 .[29] J.-P. Caltagirone, S. Vincent, On primitive formulation in fluid mechanics and fluid-structureinteraction with constant piecewise properties in velocity-potentials of acceleration, ActaMechanica 231 (6) (2020) 2155–2171. doi:10.1007/s00707-020-02630-w .[30] R. Beltman, M. Anthonissen, B. Koren, Conservative polytopal mimetic discretization ofthe incompressible Navier-Stokes equations, J. Comput. Appl. Math. 340 (2018) 443–473. doi:10.1016/j.cam.2018.02.007 .[31] R. Mittal, G. Iccarino, Immersed Boundary Methods, Annu. Rev. Fluid Mech. 37 (2005)239–261. doi:10.1146/annurev.fluid.37.061903.175743 .[32] E. Arquis, J.-P. Caltagirone, Sur les conditions hydrodynamiques au voisinage d’uneinterface milieu fluide – milieu poreux : application à la convection naturelle, C.R. Acad.Sciences, IIB 299 (1) (1984) 1–4.[33] P. Angot, C.-H. Bruneau, P. Fabrie, A penalization method to take into account obstaclesin incompressible viscous flows, Numerische Mathematik 81 (4) (1999) 497–520.[34] P. Angot, J.-P. Caltagirone, P. Fabrie, A fast vector penalty-projection methodfor incompressible non-homogeneous or multiphase Navier-Stokes problems, AppliedMathematics Letters 25 (11) (2012) 1681–1688. doi:10.1016/j.aml.2012.01.037 .[35] J.-P. Caltagirone, S. Vincent, A Kinematics Scalar Projection Method (KSP) forIncompressible Flows with Variable Density., Open Journal of Fluid Dynamics, doi:10.4236/ojfd.2015.52019 5 (2015) 171–182.[36] M. Tavares, D.-A. Koffi-Bi, E. Chenier, S. Vincent, A two-dimensional second orderconservative front-tracking method with an original marker advection approach based onjump relations, Commun. Comput. Phys. 27 (5) (2020) 1550–1589. doi:10.4208/cicp.OA-2019-0028 .[37] S. Popinet, Numerical models of surface tension, Annual Review of Fluid Mechanics 50 (1)(2018) 49–75. doi:10.1146/annurev-fluid-122316-045034 .[38] F. Denner, B. van Wachem, Numerical time-step restrictions as a result of capillary waves,Journal of Computational Physics 285 (2015) 24–40. doi:10.1016/j.jcp.2015.01.021 . 3339] S. Hysing, S. Turek, D. Kuzmin, N. Parolini, E. Burman, S. Ganesan, L. Tobiska,Quantitative benchmark computations of twodimensional bubble dynamics, Int. J. Numer.Meth. Fluids 60 (2009) 1259–1288. doi:10.1002/FLD.1934 .[40] F. Ilinca, K. R. Yu, B. Blais, The effect of viscosity on free surface flow inside an angularlyoscillating rectangular tank, Computers & Fluids 183 (2019) 160–176. doi:10.1016/j.compfluid.2019.02.021 .[41] J. Hadamard, Mouvement permanent lent d’une sphère liquide et visqueuse dans un liquidevisqueux, C. Rendus Acad. Sci. 152 (1911) 1735–1738.[42] L. Rayleigh, On the stability, or instability, of certain fluid motions, Proceedings of theLondon Mathematical Society s1-11 (1) (1879) 57–72. doi:10.1112/plms/s1-11.1.57 .[43] W. Kelvin, J. Larmor, J. Joule, Mathematical and Physical Papers, no. vol. 3 inMathematical and Physical Papers, 1890.[44] S. Chandrasekhar, The oscillations of a viscous liquid globe, Proceedings of the LondonMathematical Society s3-9 (1) (1959) 141–149. doi:10.1112/plms/s3-9.1.141 .[45] A. Prosperetti, Free oscillations of drops and bubbles: The initial-value problem, J. FluidMech. 100 (2) (1980) 333–347. doi:10.1017/S0022112080001188 .[46] A. Prosperetti, Linear oscillations of constrained drops, bubbles, and plane liquid surfaces,Physics of Fluids 24 (3) (2012) 032109. doi:10.1063/1.3697796 .[47] D. Plümacher, M. Oberlack, Y. Wang, M. Smuda, On a non-linear droplet oscillationtheory via the unified method, Physics of Fluids 32 (6) (2020) 067104. doi:10.1063/5.0007341doi:10.1063/5.0007341