Multisymplectic variational integrators for barotropic and incompressible fluid models with constraints
MMultisymplectic variational integratorsfor barotropic and incompressible fluid models withconstraints
Fran¸cois Demoures and Fran¸cois Gay-Balmaz February 23, 2021
Abstract
We present a structure preserving discretization of the fundamental spacetimegeometric structures of fluid mechanics in the Lagrangian description in 2D and 3D.Based on this, multisymplectic variational integrators are developed for barotropicand incompressible fluid models, which satisfy a discrete version of Noether theorem.We show how the geometric integrator can handle regular fluid motion in vacuumwith free boundaries and constraints such as the impact against an obstacle of afluid flowing on a surface. Our approach is applicable to a wide range of modelsincluding the Boussinesq and shallow water models, by appropriate choice of theLagrangian.
This paper presents a multisymplectic variational integrator for barotropic fluids andincompressible fluids with free boundaries in the Lagrangian description. The integratoris derived from a spacetime discretization of the Hamilton principle of fluid dynamics andis based on a discrete version of the multisymplectic geometric formulation of continuummechanics. As a consequence of its variational nature, the resulting scheme preservesexactly the momenta associated to symmetries, it is symplectic in time, and energy iswell conserved. In addition to its conservative poperties, the variational scheme can benaturally extended to handle constraints, such as the impact against an obstacle of fluidflowing on a surface, by augmenting the discrete Lagrangian with penalty terms.Multisymplectic geometry is the natural geometric setting for classical field theoriesand is the appropriate spacetime extension of the symplectic formulation of classicalmechanics. Important properties of Lagrangian and Hamiltonian systems in classicalmechanics, such as the symplecticity of the flow and the preservation of the momentum EPFL, Doc & Postdoc Alumni. Av Druey 1, Lausanne 1018, Switzerland [email protected] CNRS & ´Ecole Normale Sup´erieure, Laboratoire de M´et´eorologie Dynamique, Paris, France. [email protected] a r X i v : . [ m a t h . NA ] F e b aps associated to symmetries, have corresponding statements for field theories thatare intrinsically formulated via multisymplectic geometry. These are the multisymplecticform formula and the covariant Noether theorem for the solution of Euler-Lagrangefield equations. Of particular importance in these formulations are the Cartan forms associated to the Lagrangian density of the theory.Multisymplectic variational integrators were developed in [24] via a spacetime dis-cretization of the Hamilton principle of field theories, which results in numerical schemesthat satisfy a discrete version of the multisymplectic form formula and a discrete covari-ant Noether theorem. The discrete framework also allows the definitions of the conceptsof discrete Cartan forms and discrete covariant momentum maps . Other approaches tomultisymplectic integrators have been also developed, in, for example, [3]. We refer to[21, 8, 7, 9, 6] for the development of multisymplectic variational integrators for severalmechanical systems of interest in engineering. Examples include the simulation of thedynamics of rotor blades via asynchronous variational integrators where it is necessaryto compute accurate solutions for long periods of time, the dynamics of geometricallyexact (Cosserat) beams, or the simulation of elastodynamic frictionless impact problems.In this paper, we develop this method towards its application to compressible andincompressible fluid dynamics by using, at the continuous level, the multisymplecticvariational formulation of continuum mechanics as described in [25, 15]. The mainingredients in our discrete approach are the concepts of discrete deformation gradientand discrete Jacobian, defined both in the 2D and 3D cases. They enter in a fundamentalway in the definition of the spacetime discretized Lagrangian and they allow to exactlyimpose discrete incompressibility via an augmented Lagrangian approach. Besides itconservative properties, thanks to its variational nature, our scheme can be naturallyextended to include constraints. This is illustrated with fluid flowing or impacting on asurface.The variational discretization in this paper is carried out in the Lagrangian frame andfor fluid dynamics interpreted as a special class of field theory on spacetime. Geometricvariational discretizations for fluids have also been developed in the Eulerian descriptionand for fluid dynamics interpreted as an infinite dimensional dynamical system on dif-feomorphism groups, as opposed to the spacetime covariant description carried out here.This variational approach is based on a discretization of groups of diffeomorphisms, see[29, 30, 1, 16] for both incompressible and compressible models.This paper is a first step towards the development of dynamic mesh update from astructure preserving point of view, inspired by arbitrary Lagrangian-Eulerian methods.Several approaches have been proposed in the literature, such as [20, 11, 12, 35, 27, 13].The organization of the paper is as follows. Section § §
3. In § § § § In this section we briefly review the variational formulation of barotropic and incompress-ible fluid models in the Lagrangian (or material) description in Cartesian coordinates.This formulation is then recasted in a multisymplectic variational setting, which allowsto formulate intrinsically the Hamilton principle, the multisymplectic property of thesolutions, and the covariant Noether theorem with the help of Cartan forms. This givesthe geometric framework to be discretized in a structure preserving way later.Assume that the reference configuration of the fluid is a compact domain
B ⊂ R n with piecewise smooth boundary, and the fluid moves in the ambient space M = R n .We denote by ϕ : R × B → M the fluid configuration map, which indicates the location m = ϕ ( t, X ) at time t of the fluid particle with label X ∈ B . The deformation gradientis denoted F ( t, X ), given in coordinates by F ai = ϕ a,i , with X i , i = 1 , ..., n the Cartesiancoordinates on B and m a , a = 1 , ..., n the Cartesian coordinates on M . We assume thatthe fluid configuration is regular enough so that all the computations below are valid. A fluid is barotropic if it is compressible and the surfaces of constant pressure p andconstant density ρ coincide, i.e., we have a relation p = p ( ρ ) . (1)The internal energy W of barotropic fluids in the material description depends on thedeformation gradient F only through the Jacobian J of ϕ , given in Cartesian coordinatesby J ( t, m ) = det( F ( t, X )) , W = W ( ρ , J ), with ρ ( X ) the mass densityof the fluid in the reference configuration. The pressure in the material description is P W ( ρ , J ) = − ρ ∂W∂J ( ρ , J ) . (2)The continuity equation for mass can be written as ρ ( X ) = ρ ( t, ϕ ( t, X )) J ( t, X ) , (3)with ρ ( t, m ) the Eulerian mass density. The internal energy w ( ρ ) in the Eulerian de-scription satisfies the relation W ( ρ , J ) = w (cid:16) ρ J (cid:17) and one notes that P W = p w ◦ ϕ , with p w = ρ ∂w∂ρ the Eulerian pressure p in (1). The Lagrangian of the barotropic fluid evaluated on a fluid configuration map ϕ ( t, X )has the standard form L ( ϕ, ˙ ϕ, ∇ ϕ ) = 12 ρ | ˙ ϕ | − ρ W ( ρ , J ) − ρ Π( ϕ ) , (4)with Π a potential energy, such as the gravitational potential Π( ϕ ) = g · ϕ .Hamilton’s principle δ (cid:90) T (cid:90) B L ( ϕ, ˙ ϕ, ∇ ϕ )d t d X = 0for variations of ϕ vanishing at t = 0 , T yields the Euler-Lagrange equations ∂∂t ∂L∂ ˙ ϕ + ∂∂x i ∂L∂ϕ ,i = ∂L∂ϕ , together with the natural boundary conditions ∂L∂ϕ a,i n i δϕ a = 0 on ∂ B , for allowed variations δϕ . Here n denotes the outward pointing unit normal vector fieldto ∂ B .From the Lagrangian of the barotropic fluid (4) and the material pressure P W definedin (2) we get the barotropic fluid equations in the Lagrangian description as ρ ¨ ϕ + ∂∂x i (cid:0) P W J F − (cid:1) i = − ρ ∂ Π ∂ϕ (5)4ogether with the natural boundary conditions P W J n i ( F − ) ia δϕ a = 0 on ∂ B , (6)for all allowed variations δϕ . For instance for a free boundary problem, the variations δϕ are arbitrary on ∂ B , hence the boundary condition (6) yields the zero pressure condition P W | ∂ B = 0 . (7)Boundary conditions with surface tension can be deduced from the Hamilton principleby adding an area term in the Lagrangian, see [17].Using the relations ˙ ϕ = u ◦ ϕ , P W = p w ◦ ϕ , and ρ = ( ρ ◦ ϕ ) J , between Lagrangianand Eulerian quantities, one deduces from (5) the familiar Eulerian form of barotropicfluids as ρ ( ∂ t u + u · ∇ u ) = −∇ p w − ρ ∇ Π , ∂ t ρ + div( ρu ) = 0 . Let us consider the following general barotropic expression for the internal energy andpressure w ( ρ ) = Aγ − ρ γ − + Bρ − , p w ( ρ ) = Aρ γ − B, (8)for constants A , B , and adiabatic coefficient γ , see [5]. The material internal energy tobe used in the Lagrangian (4) is W ( ρ , J ) = Aγ − (cid:18) Jρ (cid:19) − γ + B (cid:18) Jρ (cid:19) . (9)For an isentropic perfect gas we have B = 0.In our tests, we shall use the expression (8) for the treatment of an isentropic perfectgas, where the value of the constant B (cid:54) = 0 does not affect the dynamics, while it allowsto naturally impose from (7) the boundary condition P | ∂ B = B, with P = A (cid:0) ρ J (cid:1) γ the pressure of the isentropic perfect gas. This is crucial for thediscretization, since it allows to find the appropriate discretization of the boundarycondition directly from the boundary terms of the discrete variational principle.The rotating shallow water model can also be recasted in the formulation above, inwhich case the variable ρ is interpreted as the water depth in the reference configuration.The Lagrangian is L ( ϕ, ˙ ϕ, ∇ ϕ ) = 12 ρ | ˙ ϕ | + ρ ˙ ϕ · R ( ϕ ) − ρ W ( ρ , J ) , where R is the vector potential of the angular velocity of the Earth and W is chosen as W ( ρ , J ) = g ρ J . .2 Incompressible fluid models Incompressible models are obtained by inserting the constraint J = 1 in the Hamiltonprinciple as δ (cid:90) T (cid:90) B (cid:0) L ( ϕ, ˙ ϕ, ∇ ϕ ) + λ ( J − (cid:1) d t d X = 0 , (10)where λ ( t, X ) is the Lagrange multiplier. With the Lagrangian (4), this results in thesystem ρ ¨ ϕ + ∂∂x i (cid:0) ( P W + λ ) J F − (cid:1) i = − ρ ∂ Π ∂ϕ , J = 1 . (11)With the relations ˙ ϕ = u ◦ ϕ , P W = p w ◦ ϕ , and ρ = ( ρ ◦ ϕ ) J , we get from (11) thefamiliar Eulerian formulation ρ ( ∂ t u + u · ∇ u ) = −∇ ( p w + p λ ) − ρ ∇ Π , div u = 0 , ∂ t ρ + u · ∇ ρ = 0 . In this case p w + p λ is determined from the incompressibility constraint via a Poissonequation. The
Boussinesq model is obtained from the Hamilton principle with incompressibilityconstraint (10) by interpreting ρ as the buoyancy in the reference configuration andtaking the Lagrangian L ( ϕ, ˙ ϕ, ∇ ϕ ) = 12 | ˙ ϕ | − ρ ϕ · g (12)with gravitational acceleration vector g . For the nonhomogeneous Euler fluid , the La-grangian is the kinetic energy L ( ϕ, ˙ ϕ, ∇ ϕ ) = 12 ρ | ˙ ϕ | , (13)for some non-constant density ρ ( X ). For the ideal fluid , one takes L ( ϕ, ˙ ϕ, ∇ ϕ ) = 12 | ˙ ϕ | (14)in (10), which gives ¨ ϕ + ∂∂x i (cid:0) λJ F − (cid:1) i = 0 , J = 1 , (15)and hence ∂ t u + u · ∇ u = −∇ p λ , div u = 0 is obtained in the Eulerian formulation. In this paragraph, we briefly review the geometric variational framework of classicalfield theory, as it applies to continuum mechanics, following [25]. This setting will bediscretized in a structure preserving way which allows the identification of the notion ofdiscrete multisymplecticity, discrete momentum map, and discrete Noether theorems.6 .3.1 Configuration bundle, jet bundle, and Lagrangian density
The geometric formulation of classical field theories starts with the identification of theconfiguration bundle of the theory, denoted π Y , X : Y → X , such that the fields ϕ of thetheory are sections of this fiber bundle, i.e., they are smooth maps ϕ : X → Y such that π Y , X ◦ ϕ = id X , where id X denotes the identity map on X . We assume dim X = n + 1and denote by x µ , µ = 0 , , , ..., n , the coordinates on X . The fiber coordinates on Y are y a , a = 1 , ..., N , hence coordinates on the manifold Y are ( x µ , y a ), µ = 0 , ..., n , a = 1 , ..., N . While the configuration bundle for continuum mechanics is a trivial bundle,it is advantageous to use the general setting of fiber bundles since it allows to efficentlyparticularise to continuum mechanics the intrinsic geometric formulation and structuresof field theories.The first jet bundle of the configuration bundle π X , Y : Y → X is the field theoreticanalogue of the tangent bundle of classical mechanics, i.e., its fiber at y contains thefirst derivatives ϕ a,µ ( x ) of a field ϕ at x with ϕ ( x ) = y . It is defined as the fiber bundle π Y ,J Y : J Y → Y over Y , whose fiber at y ∈ Y consists of linear maps γ : T x X → T y Y satisfying T π Y , X ◦ γ = id T x X , where x = π X , Y ( y ). The induced coordinates on the fiberof J Y → Y are denoted v aµ . We note that J Y can also be regarded as the total spaceof a bundle over X , namely π X ,J Y := π X , Y ◦ π Y ,J Y : J Y → X . Natural coordinateson the manifold J Y are hence ( x µ , y a , v aµ ), µ = 0 , ..., n , a = 1 , ..., N .The derivative of a field ϕ can be regarded as a section of π X ,J Y : J Y → X , bywriting x ∈ X (cid:55)→ j ϕ ( x ) := T x ϕ ∈ J ϕ ( x ) Y , with T x ϕ : T x X → T ϕ ( x ) Y the tangent map(or first derivative) of ϕ . The section j ϕ is called the first jet extension of ϕ and is theintrinsic object corresponding to the value of a field and of its first derivatives, at thepoints in X . In the natural coordinates ( x µ , y a , v aµ ) of J Y , the first jet extension reads j ϕ : x µ (cid:55)→ ( x µ , ϕ a ( x ) , ϕ a,µ ( x )).A Lagrangian density is a smooth bundle map L : J Y → Λ n +1 X over X , whereΛ n +1 X → X is the vector bundle of ( n + 1)-form on X . In coordinates we write L ( j ϕ ( x )) = L ( x µ , ϕ a , ϕ a,µ ) d n +1 x . The associated action functional is S ( ϕ ) := (cid:90) X L ( j ϕ ( x )) . (16) For continuum mechanics, the configuration bundle is the trivial fiber bundle Y = M × X → X , with X = R × B , where B is the reference configuration of the continuum and M is the ambient space,see the beginning of §
2. We have the equalities x = ( t, X ) and y = ( x, m ) = ( t, X, m )between the variables of the general theory and those of continuum mechanics.A section of this bundle is a map ϕ : X → X × M , whose first component is id X .It is canonically identified with a map ϕ : X = R × B → M referred to as the fluidconfiguration map above. 7he first jet bundle is canonically identified with the vector bundle L ( T X , T M ) →Y = X × M , whose fiber at y = ( x, m ) is the vector space L ( T x X , T m M ) of linear mapsfrom T x X to T m M . The first jet extension is j ϕ ( t, X ) = ( ϕ ( t, X ) , ˙ ϕ ( t, X ) , ∇ ϕ ( t, X ))and the Lagrangian density reads L ( ϕ, ˙ ϕ, ∇ ϕ ) = L ( ϕ, ˙ ϕ, ∇ ϕ ) dt ∧ d n X with L the Lagrangian of barotropic fluids given in (4). Without entering into the details, we recall that the dual jet bundle J Y (cid:63) → Y , definedas the bundle of affine maps J Y → Λ n +1 X , is endowed with a canonical ( n + 1) form Θ can and a canonical multisymplectic ( n + 2) -form Ω can = − d Θ can . These are the fieldtheoretic analogue to the canonical one-form and canonical symplectic form on the phasespace (cotangent bundle of the configuration manifold) in classical mechanics. By pullingback these canonical forms with the Legendre transform F L : J Y → J Y (cid:63) of a givenLagrangian density L : J Y → Λ n +1 X , one gets the Cartan forms Θ L and Ω L on J Y ,see [18]. These forms appear naturally in the Hamilton principle, in the multisymplecticform formula, and in the Noether theorem, as will shall explain below. All these threenotions have discrete analogues, that we shall deeply use in § § L ( j ϕ ) = ( j ϕ ) ∗ Θ L , [18], the derivative of the action functional (16) takes the intrinsicform d S ( ϕ ) · V ( ϕ ) = ddε (cid:12)(cid:12)(cid:12)(cid:12) ε =0 (cid:90) X L ( j ( φ ε ◦ ϕ ))= − (cid:90) X ( j ϕ ) ∗ i j V Ω L + (cid:90) ∂ X ( j ϕ ) ∗ i j V Θ L , (17)where φ ε is the flow of a vertical vector field V on Y , i.e., T π X , Y ◦ V = 0, and j V denotes the first jet extension of V to J Y defined as j V = ddε (cid:12)(cid:12) ε =0 j φ ε . The multisymplectic form formula is a property of the solution of the Euler-Lagrangefield equations that extends the symplectic property of the solution of the Euler-Lagrangeequations of classical mechanics. It is obtained from the identity (17), by evaluating theaction functional at a solution of the Euler-Lagrange equations and taking its derivativealong variations of solutions, see [24]. Let ϕ be a solution of the Euler-Lagrange fieldequations and V , W solutions of the first variation of the Euler-Lagrange equations at ϕ . Then V , W , ϕ satisfy the multisymplectic form formula : (cid:90) ∂U ( j ϕ ) ∗ i j V i j W Ω L = 0 , (18)for all open subset U ⊂ X with with piecewise smooth boundary.8e now recall the general statement of the covariant Noether theorem . Let a Liegroup G act on Y an assume that the action covers a diffeomorphism of X . Assumethat the Lagrangian density L is G -equivariant with respect to this action, see later in § U ⊂ X withpiecewise smooth boundary, formula (17) shows that a solution of the Euler-Lagrangefield equations satisfy the covariant Noether theorem (cid:90) ∂U ( j ϕ ) ∗ J L ( ξ ) = 0 , for all ξ ∈ g , (19)where J L ( ξ ) = i j ξ Y Θ L : J Y → g ∗ ⊗ Λ n J Y is the covariant momentum map associatedto L and ξ Y is the infinitesimal generator of the Lie group action associated to the Liealgebra element ξ ∈ g . In this section we propose a multisymplectic variational discretization of fluid mechanics,by focusing on compressible barotropic models and incompressible models. We considerfree boundary fluids, as well as fluid impacting on a surface. A main step in our con-struction is the definition of discrete deformation gradient and discrete Jacobian.
We consider the geometric setting of continuum mechanics with the configuration bundle Y = X × M → X = R × B . We assume that B is a rectangle in R and take M = R . The general discrete setting is the following. One first considers a discrete parameterspace U d and a discrete base-space configuration , which is a one-to-one map φ X d : U d → φ X d ( U d ) = X d ⊂ X whose image is the discrete spacetime X d . The discrete configuration bundle is defined as π d : Y d = X d × M → X d . The discrete fields are the sections of the discrete configurationbundle, identified with maps ϕ d : X d → M . In order to describe both the discretespacetime as well as the discrete field, one introduces the discrete configuration φ d : U d → Y , from which the discrete base-space configuration and the discrete physicaldeformation are obtained as φ X d = π d ◦ φ d and ϕ d = φ d ◦ φ − X d , see Fig. 1. This setting isparticularly well adapted to situations where the discrete spacetime is also variable, see[21, 9]. 9 d = X d × M π d (cid:15) (cid:15) U d φ d (cid:61) (cid:61) φ Xd (cid:47) (cid:47) φ X d ( U d ) = X d ⊂ X ϕ d (cid:79) (cid:79) Figure 1:
Discrete configuration and discrete configuration bundle
We consider the discrete parameter space defined by U d := { , ..., j, ..., N } × B d ,where { , ..., j, ..., N } encodes an increasing sequence of time and B d parameterizes thenodes and simplexes of the discretization of B . In this paper we restrict to the case B d = { , . . . , A } × { , . . . , B } , where A and B are the number of spatial grid points.Therefore, U d = { , . . . , N }×{ , . . . , A }×{ , . . . , B } with elements denoted ( j, a, b ) ∈ U d .The discrete parameter space determines a set of parallelepipeds, denoted (cid:28) ja,b , anddefined by the following eight pairs of indices (see Fig. 2) (cid:28) ja,b = (cid:8) ( j, a, b ) , ( j + 1 , a, b ) , ( j, a + 1 , b ) , ( j, a, b + 1) , ( j, a + 1 , b + 1) , ( j + 1 , a + 1 , b ) , ( j + 1 , a, b + 1) , ( j + 1 , a + 1 , b + 1) (cid:9) , (20) j = 0 , ..., N − a = 0 , ..., A − b = 0 , ..., B −
1. The set of all such parallelepipeds isdenoted U (cid:28) d . (a,b) (a+1,b)(a,b+1) (a+1,b+1) jj+1 Figure 2:
Discrete spacetime domain U d . As recalled above, in the continuous setting, the material internal energy function W ( ρ , J ) of the barotropic fluid depends on the deformation gradient only through itsJacobian. To define the discrete deformation gradient and the discrete Jacobian, weassume that the discrete base space configuration is of the form φ X d ( j, a, b ) = s ja,b = ( t j , z ja , z jb ) ∈ R × B , (21)see Fig. 3. The discrete field ϕ d evaluated at s ja,b is denoted ϕ ja,b := ϕ d ( s ja,b ), see Fig. 4.10 j,a,b) (j+1,a,b+1) (j,a+1,b+1)(j,a+1,b)(j,a,b+1)(j+1,a+1,b) (j+1,a+1,b+1)(j+1,a,b) ∆s a ∆s b s a+1,b+1j s a+1,b+1j+1 s a+1,bj s a+1,bj+1 s a,b+1j s a,b+1j+1 s a,bj s a,bj+1 Figure 3:
On the left : discrete coordinate system.
On the right : Nodes of the mesh with Euclideancoordinates. φ a,bj φ d (a+1,b+1)(a,b) (a+1,b)(a+1,b-1)(a,b-1)(a-1,b-1)(a-1,b) (a,b+1)(a-1,b+1)1 23 4 1 23 41 23 41 23 4 φ a,b+1j φ a+1,b+1j φ a+1,bj φ a+1,b-1j φ a-1,b+1j φ a-1,bj φ a,b-1j φ a-1,b-1j F j F j F j F j Figure 4:
Discrete field φ d = ϕ d ◦ φ X d evaluated on (cid:28) ja,b , (cid:28) ja,b , (cid:28) ja,b , (cid:28) ja,b at time t j . Given a discrete base space configuration φ X d and a discrete field ϕ d , we define thefollowing four vectors F j(cid:96) ; a,b ∈ R , (cid:96) = 1 , , , j, a, b ) ∈ U d , see Fig. 4 onthe right: F j a,b = ϕ ja +1 ,b − ϕ ja,b | s ja +1 ,b − s ja,b | and F j a,b = ϕ ja,b +1 − ϕ ja,b | s ja,b +1 − s ja,b | (22) F j a,b = ϕ ja − ,b − ϕ ja,b | s ja,b − s ja − ,b | = − F j a − ,b and F j a,b = ϕ ja,b − − ϕ ja,b | s ja,b − s ja,b − | = − F j a,b − . Based on these definitions, the discrete gradient is constructed as follows.
Definition 3.1
The discrete gradient deformations of a discrete field ϕ d at the paral-lelepiped (cid:28) ja,b are the four × matrices F (cid:96) ( (cid:28) ja,b ) , (cid:96) = 1 , , , , defined at the fournodes at time t j of (cid:28) ja,b , as follows: F ( (cid:28) ja,b ) = (cid:104) F j a,b F j a,b (cid:105) , F ( (cid:28) ja,b ) = (cid:104) F j a +1 ,b F j a +1 ,b (cid:105) , F ( (cid:28) ja,b ) = (cid:104) F j a,b +1 F j a,b +1 (cid:105) , F ( (cid:28) ja,b ) = (cid:104) F j a +1 ,b +1 F j a +1 ,b +1 (cid:105) . (23) The ordering (cid:96) = 1 to (cid:96) = 4 is respectively associated to the nodes ( j, a, b ) , ( j, a, b + 1) , ( j, a + 1 , b ) , ( j, a + 1 , b + 1) , see Fig. 4 on the left.
11t is assumed that the discrete field ϕ d is such that the determinant of the discretegradient deformations are positive. Definition 3.2
The discrete Jacobians of a discrete field ϕ d at the parallelepiped (cid:28) ja,b arethe four numbers J (cid:96) ( (cid:28) ja,b ) , (cid:96) = 1 , , , , defined at the four nodes at time t j of (cid:28) ja,b asfollows: J ( (cid:28) ja,b ) = | F j a,b × F j a,b | = det (cid:0) F ( (cid:28) ja,b ) (cid:1) ,J ( (cid:28) ja,b ) = | F j a +1 ,b × F j a +1 ,b | = det (cid:0) F ( (cid:28) ja,b ) (cid:1) J ( (cid:28) ja,b ) = | F j a,b +1 × F j a,b +1 | = det (cid:0) F ( (cid:28) ja,b ) (cid:1) J ( (cid:28) ja,b ) = | F j a +1 ,b +1 × F j a +1 ,b +1 | = det (cid:0) F ( (cid:28) ja,b ) (cid:1) . (24)As a consequence, from relations (24), the variation of the discrete Jacobian is givenby δJ (cid:96) = ∂ det( F (cid:96) ) ∂ F (cid:96) : δ F (cid:96) = J (cid:96) ( F (cid:96) ) − T : δ F (cid:96) , at each (cid:28) ja,b , which is used in the derivation of the discrete Euler-Lagrange equations. Recall that the set of all parallelepipeds in the discrete parameter space is denoted U (cid:28) d .We write X (cid:28) d := φ X d (cid:0) U (cid:28) d (cid:1) the set of all parallelepipeds in X d . The discrete version of the first jet bundle is givenby J Y d := X (cid:28) d × M × ... × M (cid:124) (cid:123)(cid:122) (cid:125) → X (cid:28) d . (25)Given a discrete field ϕ d , its first jet extension is the section of (25) defined by j ϕ d ( (cid:28) ja,b ) = (cid:0) ϕ ja,b , ϕ j +1 a,b , ϕ ja +1 ,b , ϕ j +1 a +1 ,b , ϕ ja,b +1 , ϕ j +1 a,b +1 , ϕ ja +1 ,b +1 , ϕ j +1 a +1 ,b +1 (cid:1) , (26)which associates to each parallelepiped, the values of the field at its nodes. A discreteLagrangian is a map L d : J Y d → R , see [24]. The discrete Lagrangian evaluated on a discrete field is denoted as L d (cid:0) j ϕ d ( (cid:28) ) (cid:1) . We now consider the case of the barotropic fluid. We assume for simplicity thatthe mass density ρ of the fluid in the reference configuration is a constant number.The case of a Lagrangian density with a nonconstant mass density ρ is important for12pplications to stratified flows and can be easily treated by our approach. We considera class of discrete Lagrangians associated to (4) of the form L d (cid:0) j ϕ d ( (cid:28) ) (cid:1) = vol (cid:0) (cid:28) (cid:1)(cid:16) ρ K d (cid:0) j ϕ d ( (cid:28) ) (cid:1) − ρ W d (cid:0) ρ , j ϕ d ( (cid:28) ) (cid:1) − ρ Π d (cid:0) j ϕ d ( (cid:28) ) (cid:1)(cid:17) , (27)where vol (cid:0) (cid:28) (cid:1) is the volume of the parallelepiped (cid:28) ∈ X (cid:28) d . Examples of K d , W d , Π d aregiven as follows.– The discrete kinetic energy K d : J Y d → R is defined as K d (cid:0) j ϕ d ( (cid:28) ja,b ) (cid:1) := 14 a +1 (cid:88) α = a b +1 (cid:88) β = b (cid:12)(cid:12) v jα,β (cid:12)(cid:12) , (28)with v jα,β = ( ϕ j +1 α,β − ϕ jα,β ) / ∆ t j the discrete velocity.– The discrete internal energy W d : J Y d → R is defined as W d (cid:0) ρ , j ϕ d ( (cid:28) ja,b ) (cid:1) := 14 (cid:88) (cid:96) =1 W (cid:0) ρ , J (cid:96) ( (cid:28) ja,b ) (cid:1) , (29)where W is the material internal energy of the continuous model and J (cid:96) ( (cid:28) ja,b ) arethe discrete Jacobians associated to (cid:28) ja,b at time t j .– The discrete potential energy Π d : J Y d → R is given byΠ d (cid:0) j ϕ d ( (cid:28) ja,b ) (cid:1) := 14 a +1 (cid:88) α = a b +1 (cid:88) β = b Π( ϕ jα,β ) , (30)where Π is the potential energy of the continuous model. We shall focus on thegravitation potential Π( ϕ ) = g · ϕ , with gravitational acceleration vector g , inwhich case Π d (cid:0) j ϕ d ( (cid:28) ja,b ) (cid:1) = 14 a +1 (cid:88) α = a b +1 (cid:88) β = b g · ϕ jα,β . (31)We will also consider a mid-point rule discretization later in § To simplify the exposition, we assume that the discrete base space configuration is fixedand given by φ X d ( j, a, b ) = ( j ∆ t, a ∆ s , b ∆ s ), for given ∆ t , ∆ s , ∆ s , that is, we assumethat the mesh is fixed and matches with the standard basis axis of the Euclidean space Mesh deformations can be also considered in this setting and will be explored in a future work. (cid:0) (cid:28) (cid:1) = ∆ t ∆ s ∆ s in the discrete Lagrangian(27) and the mass of each 2 D cell in φ X d ( B d ) is M = ρ ∆ s ∆ s .The discrete action functional associated to L d is obtained as S d ( ϕ d ) = (cid:88) (cid:28) ∈X (cid:28) d L d (cid:0) j ϕ d ( (cid:28) ) (cid:1) = N − (cid:88) j =0 A − (cid:88) a =0 B − (cid:88) b =0 L d (cid:0) j ϕ d ( (cid:28) ja,b ) (cid:1) . (32)In order to apply the discrete Hamilton principle, we compute the variation δS d ( ϕ d ) ofthe action sum and we get N − (cid:88) j =0 A − (cid:88) a =0 B − (cid:88) b =0 (cid:104) M (cid:16) v ja,b · δϕ j +1 a,b + v ja +1 ,b · δϕ j +1 a +1 ,b + v ja,b +1 · δϕ j +1 a,b +1 + v ja +1 ,b +1 · δϕ j +1 a +1 ,b +1 (cid:17) + A ja,b · δϕ ja,b + B ja,b · δϕ ja +1 ,b + C ja,b · δϕ ja,b +1 + D ja,b · δϕ ja +1 ,b +1 (cid:105) , where we have used the following expressions of the partial derivative of L d : D L ja,b = M v ja,b D L ja,b = M v ja,b +1 D L ja,b = M v ja +1 ,b D L ja,b = M v ja +1 ,b +1 (33)and we have introduced the following notations for the other partial derivatives D L ja,b = A ja,b D L ja,b = C ja,b D L ja,b = B ja,b D L ja,b = D ja,b , (34)whose expressions are given in Appendix A.1 for an arbitrary internal energy function W . Note that D k L ja,b is the partial derivative of L d , at (cid:28) ja,b , with respect to the k th variable, in the order listed in (26).Rearranging the expression (33) we get the discrete Euler-Lagrange equations M v j − a,b + A ja,b + B ja − ,b + C ja,b − + D ja − ,b − = 0 , (35)which correspond to variations δϕ ja,b at the interior of the domain. Variations at thespatial boundary gives the boundary conditions M v j − ,b + A j ,b + C j ,b − = 0 , M v j − a, + A ja, + B ja − , = 0 ,M v j − A,b + B jA − ,b + D jA − ,b − = 0 , M v j − a,B + C ja,B − + D ja − ,B − = 0 ,M v j − , + A j , = 0 , M v j − A, + B jA − , = 0 ,M v j − ,B + C j ,B − = 0 , M v j − A,B + D jA − ,B − = 0 , (36)14hile variations at the temporal boundary gives A a,b + B a − ,b + C a,b − + D a − ,b − = 0 , M v N − a,b = 0 ,A ,b + C ,b − = 0 , A a, + B a − , = 0 ,B A − ,b + D A − ,b − = 0 , C a,B − + D a − ,B − = 0 . (37)We assume that the variations of the discrete field at the spatial boundary are arbi-trary so we get the boundary conditions (36). This corresponds to the discrete versionof the boundary condition (6). We assume that the variations at the temporal extremityvanish, hence (37) is not imposed. In a similar way with the continuous case recalled in § L and Ω L on the first jet bundle J Y ,see § L d : J Y d → R , the discrete Cartan one-forms aredefined on the discrete first jet bundle (25) asΘ p L d = D p L d d ϕ ( p ) d , p = 1 , ..., , (38)see [24, 21, 9]. In (38) we have used the notation ϕ ( p ) d ∈ (cid:8) ϕ ja,b , ϕ j +1 a,b , ϕ ja +1 ,b , ϕ j +1 a +1 ,b , ϕ ja,b +1 , ϕ j +1 a,b +1 , ϕ ja +1 ,b +1 , ϕ j +1 a +1 ,b +1 (cid:9) . (39)For the discrete Lagrangian (27) of the barotropic fluid, using (34) and (33) weget the following expressions of the discrete Cartan one-forms evaluated on the first jetextension j ϕ d ( (cid:28) ja,b ) ∈ J Y d of discrete field ϕ d :Θ L d = A ja,b d ϕ ja,b , Θ L d = M v ja,b d ϕ j +1 a,b , Θ L d = B ja,b d ϕ ja +1 ,b , Θ L d = M v ja +1 ,b d ϕ j +1 a +1 ,b , Θ L d = C ja,b d ϕ ja,b +1 , Θ L d = M v j +1 a,b +1 d ϕ j +1 a,b +1 , Θ L d = D ja,b d ϕ ja +1 ,b +1 , Θ L d = M v ja +1 ,b +1 d ϕ j +1 a +1 ,b +1 . (40)In order to present the multisymplectic form formula and the discrete covariantNoether theorem, we shall rewrite the differential of the discrete action functional (32)in an intrinsic form using the discrete Cartan one-forms. Given a vector field V d tangent15o the discrete configuration ϕ d , we consider its first jet extension j V d which attributesto the set of nodes in (cid:28) the set of values of V d on these nodes. With this definition,for a given (cid:28) ∈ X (cid:28) d , we can write the partial derivatives of L d in terms of the discreteCartan forms as D p L d (cid:0) j ϕ d ( (cid:28) ) (cid:1) · V ( p ) d = Θ p L d (cid:0) j ϕ d ( (cid:28) ) (cid:1) · j V d = (cid:104) ( j ϕ d ) ∗ (cid:0) i j V d Θ p L d (cid:1)(cid:105) ( (cid:28) ) , (41)for p = 1 , ...,
8. In the last term we have used standard notations from differentialcalculus: the notation ( j ϕ ) ∗ for the pull-back by j ϕ d of k -forms from J Y d to X (cid:28) d andthe notation i j V d for the insertion of a vector in a k -form. With these notations, thetotal derivative of the discrete action functional (32) is d S d ( ϕ d ) · V d = (cid:88) (cid:28) ∈U (cid:28) d (cid:88) p ∈ (cid:28) (cid:104) ( j ϕ d ) ∗ (cid:0) i j V d Θ p L d (cid:1)(cid:105) ( (cid:28) ) . (42)Here p ∈ (cid:28) denotes a node p of the parallelepiped (cid:28) . Such a formula is true on anysubdomain U (cid:48) d ⊂ U d , by considering the restricted action S (cid:48) d = S d | U (cid:48) d . When restricted to a solution ϕ d of (35), the total derivative of S d reads d S d ( ϕ d ) · V d = (cid:88) (cid:28) ∈ U (cid:28) d (cid:88) p ; (cid:28) ( p ) ∈ ∂ U d (cid:104) ( j ϕ d ) ∗ (cid:0) i j V d Θ p L d (cid:1)(cid:105) ( (cid:28) ) , (43)similarly on any subdomains U (cid:48) d ⊂ U d . Note the difference with formula (42). From (43)two important results are obtained:1. The discrete multisymplectic form formula .It is obtained by taking the exterior derivative of (43), evaluating it on the firstvariations V d , W d of a solution ϕ d , and using the rules of exterior differentialcalculus, which gives dd S d ( ϕ d )( V d , W d ) = (cid:88) (cid:28) ∈ U (cid:48) (cid:28) d (cid:88) p ; (cid:28) ( p ) ∈ ∂ U (cid:48) d (cid:104) ( j ϕ d ) ∗ (cid:0) i j V d i j W d Ω p L d (cid:1)(cid:105) ( (cid:28) ) = 0 , (44)for any subdomains U (cid:48) d ⊂ U d . Here Ω p L d = − d Θ p L d , p = 1 , ..., discrete Car-tan 2-forms on J Y d , see [24]. This is the discrete version of the multisymplecticform formula (18). It extends to spacetime discretization, the symplectic propertyof variational integrators, [26]. This formula encodes a discrete version of the reci-procity theorem of continuum mechanics, as well as discrete time symplecticity ofthe solution flow, see [21]. 16. The discrete covariant Noether theorem .Consider an action Φ : G × M → M of a Lie group G on M . For ξ ∈ g , the Liealgebra of G , we denote by ξ M the infinitesimal generator of the action, i.e. thevector field on M defined by ξ M ( m ) := ddε (cid:12)(cid:12)(cid:12)(cid:12) ε =0 Φ exp( εξ ) ( m ) , for every m ∈ M . Assume that L d is G -invariant with respect this action. As aconsequence, the discrete action is also G -invariant and we get d S (cid:48) d ( ϕ d ) · ξ M ( ϕ d ) = 0 for all ξ ∈ g . (45)From (43), it follows d S (cid:48) d ( ϕ d ) · ξ M ( ϕ d ) = (cid:88) (cid:28) ∈U (cid:48) (cid:28) d (cid:88) p ; (cid:28) ( p ) ∈ ∂ U (cid:48) d (cid:104) ( j ϕ d ) ∗ (cid:10) J p L d , ξ (cid:11)(cid:105) ( (cid:28) ) = 0 , (46)for every ξ ∈ g , where the discrete covariant momentum maps are defined by J p L d : J Y d → g ∗ , (cid:104) J p L d , ξ (cid:105) := i ξ J Y d Θ p L d , ξ ∈ g , p = 1 , . . . , . (47)In (47) ξ J Y d is the infinitesimal generator of the action of G induced on J Y d bythe action Φ on M . It is given at each j ϕ d ( (cid:28) ja,b ) ∈ J Y d by ξ J Y d (cid:0) j ϕ d ( (cid:28) ja,b ) (cid:1) = (cid:16) (cid:28) ja,b , ξ M ( ϕ ja,b ) , ξ M ( ϕ j +1 a,b ) , ξ M ( ϕ ja +1 ,b ) , ξ M ( ϕ j +1 a +1 ,b ) ,ξ M ( ϕ ja,b +1 ) , ξ M ( ϕ j +1 a,b +1 ) , ξ M ( ϕ ja +1 ,b +1 ) , ξ M ( ϕ j +1 a +1 ,b +1 ) (cid:17) . From (46), we thus obtain the discrete covariant Noether theorem (cid:88) (cid:28) ∈U (cid:48) (cid:28) d (cid:88) p ; (cid:28) ( p ) ∈ ∂ U (cid:48) d J p L d ( (cid:28) ) = 0 , (48)for every subdomain U (cid:48) d ⊂ U d and for ϕ d a solution of the discrete Euler-Lagrangeequations. This is the discrete version of the covariant Noether theorem (19).We refer to [21], [8], [9] for more explanations concerning discrete conservation lawsfor multisymplectic variational discretizations. In absence of the gravitation potential, the discrete Lagrangian (27) is invariant underrotation and translation, i.e., the action of the special Euclidean group SE (2). Thisfollows from inspection of the expressions (28) and (29), and the expression of the discreteJacobian. 17rom this invariance, the discrete covariant Noether theorem (48) is satisfied withthe discrete covariant momentum maps J p L d : J Y d → se (2) ∗ given by J p L d (cid:0) j ϕ d ( (cid:28) ja,b ) (cid:1) = (cid:0) ϕ ( p ) × D p L ja,b , D p L ja,b (cid:1) , p = 1 , ..., . (49)A consequence of this discrete covariant Noether theorem is the conservation of the classical discrete momentum map given in terms of the discrete covariant momentummap as J jd = J d ( ϕ j , ϕ j +1 ) = A − (cid:88) a =0 B − (cid:88) b =0 (cid:0) J L d + J L d + J L d + J L d (cid:1) = − A − (cid:88) a =0 B − (cid:88) b =0 (cid:0) J L d + J L d + J L d + J L d (cid:1) , (50)i.e., J j +1 d = J jd . On the left hand side ϕ j = { ϕ ja,b | ≤ a ≤ A − , ≤ b ≤ B − } isthe collection of all positions at time t j . On the right hand sides each of the discretemomentum maps J p L d are evaluated on j ϕ d ( (cid:28) ja,b ). We refer to [8] for details regardingthe link between discrete classical and discrete covariant momentum maps underlyingformulas like (50). Boundary conditions play an important role in this correspondence.From (49), (50), and Appendix A.1, we get the expression J jd = A − (cid:88) a =0 B − (cid:88) b =0 J r (cid:0) j ϕ d ( (cid:28) ja,b ) (cid:1) A − (cid:88) a =0 B − (cid:88) b =0 J l (cid:0) j ϕ d ( (cid:28) ja,b ) (cid:1) with J r (cid:0) j ϕ d ( (cid:28) ja,b ) (cid:1) = a +1 (cid:88) α = a b +1 (cid:88) β = b ϕ jα,β × M v jα,β ∈ R , J l (cid:0) j ϕ d ( (cid:28) ja,b ) (cid:1) = a +1 (cid:88) α = a b +1 (cid:88) β = b M v jα,β ∈ R . (51) In this section we adapt the multisymplectic variational integrator obtained above tothe case of incompressible models.
Equality constraint.
As recalled in § J = 1. In thediscrete case, one similarly adds to the discrete action (32) the corresponding Lagrangemultiplier term to get (cid:98) S d ( ϕ d , λ d ) = (cid:88) (cid:28) ∈X (cid:28) d (cid:104) L d (cid:0) j ϕ d ( (cid:28) ) (cid:1) + (cid:88) (cid:96) =1 λ (cid:96)d ( (cid:28) ) (cid:0) J (cid:96) ( (cid:28) ) − (cid:1)(cid:105) , (52)which imposes the equality constraint J (cid:96) ( (cid:28) ) = 1 for the discrete Jacobian, for all paral-lelepiped (cid:28) and all (cid:96) = 1 , , ,
4. The critical point condition associated to (52) reads ∇ ϕ d (cid:98) S d ( ϕ d , λ d ) = 0 and J d ( (cid:28) ) − (cid:28) ∈ X (cid:28) d . (53)18t is well-known, [33, p.187], that if ϕ d is a local optimal solution of the function S d ( ϕ d ) in (32) restricted to the discrete incompressibility equality constraint C d = { ϕ d | J (cid:96) ( (cid:28) ) = 1 , ∀ (cid:28) ∈ X (cid:28) d } , then there must be a Lagrangian multiplier λ d such that (53)holds. However, solving the equations in (53) may not be practical, because inequalityconstraints also naturally appear, as we will see in the following examples. Let us thusconsider a constraint set C d associated to the equality constraint J (cid:96) = 1 and to inequalityconstraints g i ≤ i = 1 , ...m , i.e., C d = { ϕ d | J (cid:96) ( (cid:28) ) = 1 and g i ( (cid:28) ) ≤ i = 1 , ...m } . (54)Under appropriate conditions, see [34], generalizations of the Lagrangian multiplier ruleallow to find the critical points of the action S d defined in (32) under constraints of theform (54), see [9] for an application to variational integrators. In particular, we have thefollowing necessary condition for ϕ d ∈ C d to be locally optimal: −∇ S d ( ϕ d ) ∈ N C d ( ϕ d ),where N C d ( ϕ d ) is the normal cone to C d at ϕ d , which can be viewed as a special case ofthe calculus of subgradients, i.e., N C d ( ϕ d ) = ∂I C d ( ϕ d ), where I C d is the indicator functionof C d .When C d and S d are convex, the locally optimal condition is sufficient for ϕ d to beglobally optimal. If the ambient space is M = R n , this relation reduces to ∇ S d ( ϕ d ) + λ (cid:96) ∇ J (cid:96) ( ϕ d ) + λ ∇ g ( ϕ d ) + .... + λ m ∇ g m ( ϕ d ) = 0 , (55)for λ (cid:96) ∈ R and where λ i ≥
0, for i = 1 , ..., m , is non-vanishing only when g i ( ϕ d ) = 0.When S d is not convex it is difficult, on a practical viewpoint, to find the global optimalamong the set of local optimals, see e.g., [9].Also, given the examples that we will study, where C d is convex, instead of solv-ing our problem with the Lagrangian multiplier approach we will introduce quadraticpenalty functions rα ( ϕ d ) associated with (54), with penalty parameter r , which may beconsidered as an approximation of the indicator function I C d , see [28, p.280]. Moreover,if we suppose that for each r there exists a solution ϕ r ∈ M to the problem to minimize S d ( ϕ d ) + r α ( ϕ d ) with ϕ d ∈ M , and that the sequence { ϕ r } is contained in a compactsubset of M , we know (see [2, p.477]) that the limit ϕ d of any convergent subsequence of { ϕ r } when r → ∞ is an optimal solution to the original problem. If C d is nonconvex, alarge enough penalty parameter r must be used to get sufficiently close to an optimal so-lution. In this case computational difficulties could appear to solve the penalty problemand an augmented Lagrangian penalty function can be considered, which enjoys severaladvantegous properties, see, e.g., [33, 34, 2]. Penalty method.
Given the discrete action S d defined in (32), the Hamilton principlesubject to constraints is approximated by a penalty scheme where one seeks the criticalpoints of the action (cid:101) S d ( ϕ d ) = S d ( ϕ d ) − (cid:88) (cid:28) ∈X (cid:28) d vol( (cid:28) ) (cid:0) Φ d (cid:0) j ϕ d ( (cid:28) ) (cid:1) + Φ d (cid:0) j ϕ d ( (cid:28) ) (cid:1) + ... + Φ dm (cid:0) j ϕ d ( (cid:28) ) (cid:1)(cid:1) , (56)19ith quadratic penalty term associated to the incompressibility (equality) constraintΦ d (cid:0) j ϕ d ( (cid:28) ) (cid:1) := 14 (cid:88) (cid:96) =1 r (cid:0) J (cid:96) ( (cid:28) ) − (cid:1) , (57)where r is the penalty parameter, and with quadratic penalty terms Φ di (cid:0) j ϕ d ( (cid:28) ) (cid:1) , i = 2 , ...m , associated to inequality constraints. To test the convergence of our multisymplectic integrator, we will use an implicit inte-grator obtained from the Lagrangian (27) discretized through the mid-point rule, thatis, the discrete internal energy W d (cid:0) ρ , j ϕ d ( (cid:28) ja,b ) (cid:1) is evaluated at time ( t j + t j +1 ) /
2. Inthis case the four vectors (22) at each node ( j, a, b ) ∈ U d are now defined by F j +1 / a,b = ( ϕ ja +1 ,b + ϕ j +1 a +1 ,b ) − ( ϕ ja,b + ϕ j +1 a,b )2∆ s , F j +1 / a,b = ( ϕ ja,b +1 + ϕ j +1 a,b +1 ) − ( ϕ ja,b + ϕ j +1 a,b )2∆ s F j +1 / a,b = ( ϕ ja − ,b + ϕ j +1 a − ,b ) − ( ϕ ja,b + ϕ j +1 a,b )2∆ s , F j +1 / a,b = ( ϕ ja,b − + ϕ j +1 a,b − ) − ( ϕ ja,b + ϕ j +1 a,b )2∆ s . (58) The variation δS d ( ϕ d ) of the action sum is found as N − (cid:88) j =0 A − (cid:88) a =0 B − (cid:88) b =0 (cid:20)(cid:18) M v ja,b + A ja,b (cid:19) · δϕ j +1 a,b + (cid:18) M v ja +1 ,b + B ja,b (cid:19) · δϕ j +1 a +1 ,b + (cid:18) M v ja,b +1 + C ja,b (cid:19) · δϕ j +1 a,b +1 + (cid:18) M v ja +1 ,b +1 + D ja,b (cid:19) · δϕ j +1 a +1 ,b +1 + (cid:18) − M v ja,b + A ja,b (cid:19) · δϕ ja,b + (cid:18) − M v ja +1 ,b + B ja,b (cid:19) · δϕ ja +1 ,b + (cid:18) − M v ja,b +1 + C ja,b (cid:19) · δϕ ja,b +1 + (cid:18) − M v ja +1 ,b +1 + D ja,b (cid:19) · δϕ ja +1 ,b +1 (cid:21) , with coefficients A ja,b , B ja,b , C ja,b , D ja,b given in Appendix A.1 for an arbitrary internalenergy function W . It yields the implicit discrete Euler-Lagrange equations Mv j − a,b − Mv ja,b + A ja,b + A j − a,b + B ja − ,b + B j − a − ,b + C ja,b − + C j − a,b − + D ja − ,b − + D j − a − ,b − = 0 , (59) with the spatial boundary conditions M v j − ,b − M v j ,b + A j ,b + A j − ,b + C j ,b − + C j − ,b − = 0 ,M v j − a, − M v ja, + A ja, + A j − a, + B ja − , + B j − a − , = 0 ,M v j − A,b − M v jA,b + B jA − ,b + B j − A − ,b + D jA − ,b − + D j − A − ,b − = 0 ,M v j − a,B − M v ja,B + C ja,B − + C j − a,B − + D ja − ,B − + D j − a − ,B − = 0 ,M v j − , − M v j , + A j , + A j − , = 0 , M v j − A, − M v jA, + B jA − , + B j − A − , = 0 ,M v j − ,B − M v j ,B + C j ,B − + C j − ,B − = 0 , M v j − A,B − M v jA,B + D jA − ,B − + D j − A − ,B − = 0 . (60) The corresponding temporal boundary conditions can be computed similarly.20 .3 Numerical simulation
We evaluate the properties of the proposed multisymplectic integrator for barotropicand incompressible ideal fluid models with the case of a free boundary fluid and withthe case of a fluid flowing on a surface and impacting an obstacle. In the two cases wepresent an explicit integrator, while we consider an implicit integrator (mid-point rule)for the convergence tests.
Consider a barotropic fluid with properties ρ = 997 kg / m , γ = 6, A = ˜ Aρ − γ with˜ A = 3 . × Pa, and B = 3 . × Pa. The size of the discrete referenceconfiguration at time t is 1m × s = ∆ s = 0 . r = 0) and the incompressible casewith penalty parameters r = 10 and r = 10 . The time-steps are ∆ t = 10 − when r ∈ { , } , and ∆ t = 5 × − when r = 10 . Initial perturbations (tiny compression)are applied at time t on nodes (4 ,
0) and (5 , W plays no role since its effect is absorbed into the gradient of the pressure. In the discretecase, when using the penalty method for incompressible fluids it is advantageous toinclude the internal energy of the isentropic perfect fluid, as the case W = 0 needs todeal with a much higher penalty term.Figure 5: Left to right : barotropic and incompressible ideal fluid models ( r = 10 and r = 10 ). Topto bottom : after 0 .
1s and 6s.
21e observe that the compressible model exhibits enhanced deformation. The differ-ent behavior of the compressible and incompressible model will be even more noticeablein the following test, see Fig. 8 with respect to Fig. 7.The discrete Lagrangian is invariant under rotation and translation, hence from thediscrete Noether theorem the angular and linear momentum map (51) are preserved.Energy and momentum preservation is illustrated in Fig. 6.Figure 6:
Left to right : barotropic and incompressible ideal fluid models ( r = 10 and r = 10 ). Topto bottom : Relative energy and momentum map evolution during 6s.
Let us consider a fluid subject to gravity and flowing withoutfriction on a surface until it comes into contact with an obstacle. The gravitationalpotential is described by (31) with g = g E . We consider both a barotropic fluid and anincompressible ideal fluid.We impose the following constraints on the configuration:– The fluid is bounded below by a rigid surface, defined by the inequality constraintΨ ( ϕ ja,b ) ≤ ϕ ja,b .– There is a second inequality constraint Ψ ( ϕ ja,b ) ≤
0, verified for all ϕ ja,b , whichforces the fluid to stay outside of the obstacle.For barotropic fluid and incompressible ideal fluid the problems to solve are respec-tively described as follows– ( P ) Find the critical points of the action S d defined in (32) subject to the in-equality constraints Ψ α ( ϕ ja,b ) ≤ , α = 1 , , for all nodes.
22 ( P ) Find the critical points of the action S d defined in (32) subject to the equalityconstraint J (cid:96) ( (cid:28) ) = 1 , see (52) , and the inequality constraints Ψ α ( ϕ ja,b ) ≤ , α =1 , , for all nodes. We solve the previous problems via the penalty method. For problem P , we mustfind the critical points of the action (cid:101) S d ( ϕ d ) = S d ( ϕ d ) − N − (cid:88) j =0 A − (cid:88) a =0 B − (cid:88) b =0 ∆ t ∆ s ∆ s (cid:16) Φ d ( ϕ ja,b ) + Φ d ( ϕ ja,b ) (cid:17) , (61)Φ dα ( ϕ ja,b ) = 12 K α | Ψ α ( ϕ ja,b ) | with (cid:40) K α ∈ ]0 , ∞ [ if Ψ α ( ϕ ja,b ) ≥ K α = 0 if Ψ α ( ϕ ja,b ) < . (62)For problem P we add the penalty function (57), associated to the equality constraint J (cid:96) ( (cid:28) ) = 1, into the discrete action (61). Test.
Consider a barotropic fluid model with properties ρ = 997 kg / m , γ = 6, A =˜ Aρ − γ with ˜ A = 3 . × Pa, and B = 3 . × Pa. The size of the discretereference configuration at time t is 2 m × . t = 10 − and space-steps∆ s = 0 . s = 0 . K = 4 . × , K = 4 . × . For the incompressible case we considerthe penalty parameter r = 5 × .The initial motion of the fluid is only due to the gravity. There are no other pertur-bations so that there is no expansion or compression imposed in the initial conditions.The evolution in the barotropic and incompressible cases are illustrated in Fig. 7 andFig. 8, with the incompressibility conditions imposed by the penalty term.23igure 7: Barotropic fluid with contact.
Top to bottom : after 0 . . . Figure 8:
Incompressible ideal fluid with contact.
Top to bottom : after 0 . . . Left to right
Barotropic fluid and incompressible ideal fluid with contact.
Top to bottom :Relative energy and momentum map evolution during 2s.
Consider a barotropic fluid model with properties ρ = 997 kg / m , γ = 6, A = ˜ Aρ − γ with ˜ A = 3 . × Pa, and B = 3 . × Pa. The size of the discrete referenceconfiguration at time t is 0 . × . implicit integrator to studythe convergence with respect to ∆ t and ∆ s = ∆ s . Barotropic fluid motion in vacuum with free boundaries.
Given a fixed mesh, withvalues ∆ s = ∆ s = 0 . V a, = (0 , . × a ) T on theboundary b = 0 for all a ∈ { , ..., A } , and vary the time-steps as ∆ t ∈ { × − , . × − , . × − , . × − } . We compute the L -errors in the position ϕ d at time t N = 0 . ϕ d with an “exact solution” obtained with the time-step∆ t ref = 3 . × − s . That is, for each value of ∆ t we calculate (cid:107) ϕ d − ϕ ref (cid:107) L = (cid:32)(cid:88) a (cid:88) b (cid:107) ϕ Na,b − ϕ N ref; a,b (cid:107) (cid:33) / . (63)This yields the following convergence with respect to ∆ t t × − . × − . × − . × − (cid:107) ϕ d − ϕ ref (cid:107) L . × − × − . × − × − rate 1.106 0.964 1.815Given a fixed time-step ∆ t = 2 × − we impose an initial speed V a, = (0 , . × a ) T , on the boundary b = 0 for all a ∈ { , ..., A } , and vary the space-steps as ∆ s = ∆ s ∈ { . , . , . , . } . The “exact solution” is chosen with ∆ s = ∆ s = 0 . L -errors in the position ϕ d at time t N = 0 . s = ∆ s ∆ s = ∆ s . . . . (cid:107) ϕ d − ϕ ref (cid:107) L × − . × − . × − . × − rate 0.475 0.493 1.12 Impact against an obstacle of a fluid flowing on a surface.
The values of the impene-trability penalty parameters are K = K = 3 × . Given a fixed mesh, with values∆ s = ∆ s = 0 . t ∈ { × − , . × − , . × − , . × − } . Then, we compute the L -errors in the position ϕ d at time t N = 0 . ϕ d with an “exact solution”obtained with the time-step ∆ t ref = 3 . × − . We get the following convergence withrespect to ∆ t ∆ t × − . × − . × − . × − (cid:107) ϕ d − ϕ ref (cid:107) L . × − . × − . × − × − rate 1.101 1.226 1.586Similarly, given a fixed time-step ∆ t = 2 × − we repeat the same experimentwith varying values of ∆ s = ∆ s ∈ { . , . , . , . } . The “exact solution” is chosenwith ∆ s = ∆ s = 0 . L -errors in the position ϕ d at time t N = 0 . s = ∆ s ∆ s = ∆ s . . . . (cid:107) ϕ d − ϕ ref (cid:107) L . × − . × − . × − . × − rate 0.733 1.132 1.366An illustration of the test used for the numerical convergence is given in Fig. 10. Note that, we need to take care of the initial sum of momentum (cid:80) a m a V a, which must be of thesame value regardless of the number of nodes in the mesh. Barotropic fluid.
From top to bottom : motion in vacuum with free boundaries (∆ s = ∆ s =0 .
05, ∆ t = 2 × − ), and impact against an obstacle of a fluid flowing on a surface (∆ s = ∆ s = 0 . t = 2 × − ). From left to right : after 0 . . . Figure 11:
Relative error in the energy ( ET j − ET ) /ET , with ET j the total energy at time t j . Fromleft to right : motion in vacuum with free boundaries (∆ s = ∆ s = 0 .
05, ∆ t = 2 × − ), and impactagainst an obstacle of a fluid flowing on a surface (∆ s = ∆ s = 0 .
05, ∆ t = 2 × − ). In this section we indicate how the developments made in § §
3. We assume that B is a parallelepiped in R and take M = R . The discrete parameter space U d is now decomposed in a set of elements (cid:28) ja,b,c definedby 16 pairs of indices, see Fig. 12 for the eight pairs of indices in (cid:28) ja,b,c at time t j .As in (21), we consider discrete base-space configurations of the form φ X d : U d (cid:51) ( j, a, b, c ) (cid:55)→ s ja,b,c = ( t j , z ja , z jb , z jc ) ∈ X d ⊂ R × B . (64)The discrete field ϕ d evaluated at s ja,b,c is denoted ϕ ja,b,c := ϕ d ( s ja,b,c ) ∈ R .27 (a+1,b-1)(a,b-1)(a-1,b-1) φ a+1,b-1j φ a,b-1j φ a-1,b-1j (a,b,c) (a+1,b,c)1 234 57 8 (a+1,b+1,c)(a+1,b+1,c+1)(a+1,b,c+1)(a,b+1,c+1)(a,b+1,c)(a,b,c+1) φ d φ a,b,cj φ a,b,c+1j φ a,b+1,c+1j φ a,b+1,cj φ a+1,b+1,c+1j φ a+1,b+1,cj φ a+1,b,c+1j φ a+1,b,cj F j F j F j F j F j F j Figure 12:
Discrete field φ d = ϕ d ◦ φ X d evaluated on (cid:28) ja,b,c at time t j . Given a discrete base space configuration φ X d of the form (64) and a discrete field ϕ d ,we define the following six vectors F j(cid:96) ; a,b,c ∈ R , (cid:96) = 1 , ..., j, a, b, c ) ∈ U d ,see Fig. 12 on the right: F j a,b,c = ϕ ja +1 ,b,c − ϕ ja,b,c | s a +1 ,b,c − s a,b,c | , F j a,b,c = ϕ ja,b +1 ,c − ϕ ja,b,c | s a,b +1 ,c − s a,b,c | , F j a,b,c = ϕ ja,b,c +1 − ϕ ja,b,c | s a,b,c +1 − s a,b,c | , F j a,b,c = ϕ ja − ,b,c − ϕ ja,b,c | s a,b,c − s a − ,b,c | = − F j ,a − ,b,c F j a,b,c = ϕ ja,b − ,c − ϕ ja,b,c | s a,b,c − s a,b − ,c | = − F j ,a,b − ,c , F j a,b,c = ϕ ja,b,c − − ϕ ja,b,c | s a,b,c − s a,b,c − | = − F j a,b,c − . Based on these definitions, the discrete gradient is constructed as follows.
Definition 4.1
The discrete gradient deformations of a discrete field ϕ d at the element (cid:28) ja,b,c are the × matrices F (cid:96) ( (cid:28) ja,b,c ) , (cid:96) = 1 , ..., , defined at the eight nodes at time t j of (cid:28) ja,b,c as follows: F ( (cid:28) ja,b,c ) = (cid:104) F j a,b,c F j a,b,c F j a,b,c (cid:105) , F ( (cid:28) ja,b,c ) = (cid:104) F j a +1 ,b,c F j a +1 ,b,c F j a +1 ,b,c (cid:105) , F ( (cid:28) ja,b,c ) = (cid:104) F j a,b +1 ,c F j a,b +1 ,c F j a,b +1 ,c (cid:105) , F ( (cid:28) ja,b,c ) = (cid:104) F j a,b,c +1 F j a,b,c +1 F j a,b,c +1 (cid:105) , F ( (cid:28) ja,b,c ) = (cid:104) F j a +1 ,b +1 ,c F j a +1 ,b +1 ,c F j a +1 ,b +1 ,c (cid:105) , F ( (cid:28) ja,b,c ) = (cid:104) F j a,b +1 ,c +1 F j a,b +1 ,c +1 F j a,b +1 ,c +1 (cid:105) , F ( (cid:28) ja,b,c ) = (cid:104) F j a +1 ,b,c +1 F j a +1 ,b,c +1 F j a +1 ,b,c +1 (cid:105) , F ( (cid:28) ja,b,c ) = (cid:104) F j a +1 ,b +1 ,c +1 F j a +1 ,b +1 ,c +1 F j a +1 ,b,c +1 (cid:105) . (65)28 he ordering (cid:96) = 1 to (cid:96) = 8 is respectively associated to the nodes ( j, a, b, c ) , ( j, a +1 , b, c ) , ( j, a, b + 1 , c ) , ( j, a, b, c + 1) , ( j, a + 1 , b + 1 , c ) , ( j, a, b + 1 , c + 1) , ( j, a + 1 , b, c + 1) , ( j, a + 1 , b + 1 , c + 1) , see Fig. 12 on the left. Then we define the Jacobian in each node, as follows
Definition 4.2
The discrete Jacobians of a discrete field ϕ d at the element (cid:28) ja,b,c are thenumbers J (cid:96) ( (cid:28) ja,b ) , (cid:96) = 1 , ..., , defined at the eight nodes at time t j of (cid:28) ja,b,c as follows: J ( (cid:28) ja,b,c ) = ( F j a,b,c × F j a,b,c ) · F j a,b,c = det (cid:0) F ( (cid:28) ja,b,c ) (cid:1) . (66)See in § A.3 for the others Jacobian on (cid:28) ja,b,c . We can now establish the link betweenthe discrete Jacobian and the discrete gradient deformation.In terms of the discrete field ϕ d , the discrete Jacobians are J ( (cid:28) ja,b,c ) = (( ϕ ja +1 ,b,c − ϕ ja,b,c ) × ( ϕ ja,b +1 ,c − ϕ ja,b,c )) · ( ϕ ja,b,c +1 − ϕ ja,b,c ) | s a +1 ,b,c − s a,b,c || s a,b +1 ,c − s a,b,c || s a,b,c +1 − s a,b,c | J ( (cid:28) ja,b,c ) = (( ϕ ja +1 ,b +1 ,c − ϕ ja +1 ,b,c ) × ( ϕ ja,b,c − ϕ ja +1 ,b,c )) · ( ϕ ja +1 ,b,c +1 − ϕ ja +1 ,b,c ) | s a +1 ,b +1 ,c − s a +1 ,b,c || s a,b,c − s a +1 ,b,c || s a +1 ,b,c +1 − s a +1 ,b,c | similarly for the other ones. The discrete Lagrangian for 3D barotropic fluid models has the same general form as(27), with the obvious 3D extension of formulas (28)–(31).
The discrete action functional takes the form S d ( ϕ d ) = N − (cid:88) j =0 A − (cid:88) a =0 B − (cid:88) b =0 C − (cid:88) b =0 L (cid:0) j ϕ d ( (cid:28) ja,b,c ) (cid:1) (67)and yield the discrete Euler-Lagrange equations M v ja,b,c + A ja,b,c + B ja − ,b,c + C ja,b − ,c + D ja,b,c − + E ja − ,b − ,c + F ja,b − ,c − + G ja − ,b,c − + H ja − ,b − ,c − = 0 , (68)where we have used notations analogous to (33) and (34) for the partial derivative of L d .We refer to Appendix A.4 for the expressions of A ja,b,c , ..., H ja,b,c . Boundary conditionsare deduced the discrete Hamilton principle in a similar way as it was done in (36) and(37) for the 2D case. 29 .1.5 Discrete multisymplectic form formula and discrete Noether theorem Following the general definition (38), the discrete Cartan forms evaluated at the first jetextension j ϕ d ( (cid:28) ja,b,c ) of a discrete field ϕ d are Θ L d = A ja,b,c d ϕ ja,b,c , Θ L d = M v ja,b,c d ϕ j +1 a,b,c , Θ L d = B ja,b,c d ϕ ja +1 ,b,c , Θ L d = M v ja +1 ,b,c d ϕ j +1 a +1 ,b,c , Θ L d = C ja,b,c d ϕ ja,b +1 ,c , Θ L d = M v j +1 a,b +1 ,c d ϕ j +1 a,b +1 ,c , Θ L d = D ja,b,c d ϕ ja,b,c +1 , Θ L d = M v ja,b,c +1 d ϕ j +1 a,b,c +1 , Θ L d = E ja,b,c d ϕ ja +1 ,b +1 ,c , Θ L d = M v ja +1 ,b +1 ,c d ϕ j +1 a +1 ,b +1 ,c , Θ L d = F ja,b,c d ϕ ja,b +1 ,c +1 , Θ L d = M v ja,b +1 ,c +1 d ϕ j +1 a,b +1 ,c +1 , Θ L d = G ja,b,c d ϕ ja +1 ,b,c +1 , Θ L d = M v j +1 a +1 ,b,c +1 d ϕ j +1 a +1 ,b,c +1 , Θ L d = H ja,b,c d ϕ ja +1 ,b +1 ,c +1 , Θ L d = M v ja +1 ,b +1 ,c +1 d ϕ j +1 a +1 ,b +1 ,c +1 . (69) With these forms, the discrete multisymplectic form formula and conservation lawsin the presence of a symmetry group (discrete Noether theorem) can be derived in asimilar way as it was done in 2D in § Exactly as in § SE (3) invariant and hence the discretecovariant Noether theorem holds with the covariant discrete momentum maps J p L d : J Y d → se (3) ∗ , p = 1 , ...,
16. From this, the discrete momentum map J jd = J d ( ϕ j , ϕ j +1 ) = A − (cid:88) a =0 B − (cid:88) b =0 C − (cid:88) c =0 (cid:0) J L d + J L d + J L d + J L d + J L d + J L d + J L d + J L d (cid:1) = − A − (cid:88) a =0 B − (cid:88) b =0 C − (cid:88) c =0 (cid:0) J L d + J L d + J L d + J L d + J L d + J L d + J L d + J L d (cid:1) , (70) is preserved, as explained in § J jd = A − (cid:88) a =0 B − (cid:88) b =0 C − (cid:88) a =0 J r (cid:0) j ϕ d ( (cid:28) ja,b ) (cid:1) A − (cid:88) a =0 B − (cid:88) b =0 C − (cid:88) a =0 J l (cid:0) j ϕ d ( (cid:28) ja,b ) (cid:1) J r (cid:0) j ϕ d ( (cid:28) ja,b,c ) (cid:1) = a +1 (cid:88) α = a b +1 (cid:88) β = b c +1 (cid:88) α = c ϕ jα,β,γ × (cid:18) M v jα,β,γ (cid:19) ∈ R J l (cid:0) j ϕ d ( (cid:28) ja,b,c ) (cid:1) = a +1 (cid:88) α = a b +1 (cid:88) β = b c +1 (cid:88) α = c M v jα,β,γ ∈ R . As done in 2D (see § J = 1 we consider apenalty function Φ d (cid:0) j ϕ d ( (cid:28) ) (cid:1) := 18 (cid:88) (cid:96) =1 r (cid:0) J (cid:96) ( (cid:28) ) − (cid:1) , (71)where r is the penalty parameter. In this section we illustrate the performance of an explicit-in-time integrator in 3D , asit was done in § Consider a barotropic fluid with properties ρ = 997 kg / m , γ = 6, A = ˜ Aρ − γ with˜ A = 3 . × Pa, and B = 3 . × Pa. The size of the mesh at time t is2 m × × s = ∆ s = ∆ s = 0 .
333 m. We consider both the compressiblebarotropic fluid and the incompressible case with penalty parameter r = 10 and r = 10 .The time step is ∆ t = 10 − . There are no exterior forces.Initial perturbations are applied at time t on nodes (1 , ,
1) and (1 , , Note that with the multisymplectic variational integrators we can equally move in time and in space,see [8].
Left to right : Discrete barotropic and incompressible ideal fluid model ( r = 10 and r = 10 ). Top to bottom : after 0 .
1s and 2s.
Figure 14:
Left to right : barotropic and incompressible ideal fluid model ( r = 10 and r = 10 ). Topto bottom : Relative energy and momentum map evolution during 3s.
The main interest of this test in vacuum with free boundaries is to exhibit the perfectpreservation of the symmetries, see the figures above.
As explain in § P ), resp., ( P ) denote the problem to solve for barotropic fluid, resp., incom-pressible ideal fluid. These two problems are already described in § ρ = 997 kg / m , γ = 7, A = ˜ Aρ − γ with˜ A = 3 . × Pa, and B = 3 . × Pa. The size of the discrete referenceconfiguration at time t is 1 . × × . t = 5 × − and space-steps ∆ s = 0 . s = 0 . s = 0 . K = K = 5 × . We consider both the compressible barotropic fluidand the incompressible case with penalty given by r = 10 and r = 10 .As in 2D, the initial motion of the fluid is only due to the gravity.Figure 15: Barotropic fluid model with impact after 0 . . . . Fluid impact.
Left to right : From less ( r = 10 ) to more ( r = 10 ) incompressibility. Top tobottom : after 0 . . . . Note the large differences in behavior between the three tests that correspond tobarotropic fluid (Fig. 15) and to fluids which are more or less incompressible (Fig. 16).34 steps -300-200-1000100200300 J J J J J J steps -300-200-1000100200300 J J J J J J steps -300-200-1000100200300 J J J J J J Figure 17:
From left to right : Barotropic, incompressible ideal ( r = 10 ), incompressible ideal ( r = 10 )fluid impact model. From top to bottom : Relative energy and momentum map evolution.
We observe that three components of the momentum map are preserved until thecontact with the obstacle. These three components correspond to the symmetries ofthe gravitational term, given by the three dimensional subgroup of SE (3) consisting ofrotations around the vertical axis and translations in the horizontal plane. After contactonly one symmetry with respect to one axis of translation is conserved, namely, thetranslation parallel to the obstacle wall. After impact, the total energy becomes moreunstable, as observed in the 2D test, see Fig. 9.Results concerning impact between fluid and solid are promising, however furtherinvestigations in numerics are necessary. In particular, we must develop an implicitintegrator in order to increase the time-step and the performance of the integrator whenthe stability of the flow begin to be lost (e.g., when the flow is perturbed by the impact). Consider a barotropic fluid model with properties ρ = 997 kg / m , γ = 6, A = ˜ Aρ − γ with ˜ A = 3 . × Pa, and B = 3 . × Pa. The size of the discrete referenceconfiguration at time t is 0 . × . × . explicit integrator tostudy the convergence with respect to ∆ t and ∆ s i , i = 1 , , Barotropic fluid flowing freely over a surface.
Given a fixed mesh, with values ∆ s =∆ s = ∆ s = 0 . t ∈ { . × − , . × − , . × − , . × − } . We computethe L -errors in the position ϕ d at time t N = 0 . ϕ d with an “exactsolution” obtained with the time-step ∆ t ref = 7 . × − s . That is, for each value of ∆ t
35e calculate (cid:107) ϕ d − ϕ ref (cid:107) L = (cid:32)(cid:88) a (cid:88) b (cid:88) c (cid:107) ϕ Na,b,c − ϕ N ref; a,b,c (cid:107) (cid:33) / . (72)This yields the following convergence with respect to ∆ t ∆ t . × − . × − . × − . × − (cid:107) ϕ d − ϕ ref (cid:107) L × − . × − . × − . × − rate 1.03 1.07 1.11Given a fixed time-step ∆ t = 3 . × − , we vary the space-steps as ∆ s = ∆ s =∆ s ∈ { . , . , . , . } . The “exact solution” is chosen with ∆ s = ∆ s =∆ s = 0 . L -errors in the position ϕ d at time t N = 0 . s = ∆ s = ∆ s ∆ s = ∆ s = ∆ s . . .
05 0 . (cid:107) ϕ d − ϕ ref (cid:107) L . . . . Barotropic fluid flowing freely over a surface with ∆ t = 3 . × − . From left to right :after 0 .
25s and 0 . From top to bottom : with ∆ s i = 0 .
1m and ∆ s i = 0 . Relative error in the energy ( ET j − ET ) /ET . From left to right : (∆ s i = 0 . , ∆ t =2 . × − ), (∆ s i = 0 . , ∆ t = 3 . × − ), (∆ t = 3 . × − , ∆ s i = 0 . This paper has presented new Lagrangian schemes for the regular motion of barotropicand incompressible fluid models, which preserve the momenta associated to symmetries,up to machine precision, and satisfy the nearly constant energy property of symplecticintegrators, see Fig. 6 and 14. The schemes are derived by discretization of the geometricand variational structures underlying the spacetime formulation of continuum mechanics,seen as a particular instance of field theory. We have illustrated how this approachcan naturally accommodate incompressibility and fluid impact against an obstacle byappropriately augmenting the discrete Lagrangian and thanks to our definition of thediscrete Jacobian.An important task for the future is the study of contacts, and of friction with heatexchange, between liquid and gas, which are important phenomena encountered forinstance in ocean-atmosphere coupling.Thanks to the clear link between the expressions of the discrete Jacobian and discretegradient deformation in § § A Appendix
A.1 Derivatives of the discrete Lagrangian for 2D barotropic fluid
Explicit integrator:
The partial derivatives of the discrete Lagrangian for a 2D barotropicfluid with internal energy W ( ρ , J ) are listed below, where the discrete pressure at time t j and spatial positon (cid:96) , with the ordering (cid:96) = 1 to (cid:96) = 4 respectively associated to thenodes ( j, a, b ), ( j, a + 1 , b ), ( j, a, b + 1), ( j, a + 1 , b + 1), is defined as P (cid:96) ( (cid:28) ja,b ) = − ρ ∂W ( (cid:28) ja,b ) ∂J (cid:96) ( (cid:28) ja,b ) = P W (cid:0) ρ , J (cid:96) ( (cid:28) ja,b ) (cid:1) . P (cid:96) ( (cid:28) ja,b ) = A (cid:16) ρ J (cid:96) ( (cid:28) ja,b ) (cid:17) γ − B. (73) A ja,b = − M v ja,b + ∆ t P ( (cid:28) ja,b ) (cid:16) ϕ ja +1 ,b − ϕ ja,b +1 (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ϕ ja +1 ,b − ϕ ja +1 ,b +1 (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ϕ ja +1 ,b +1 − ϕ ja,b +1 (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | B ja,b = − M v ja +1 ,b + ∆ t P ( (cid:28) ja,b ) (cid:16) ϕ ja,b +1 − ϕ ja,b (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ϕ ja +1 ,b +1 − ϕ ja,b (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ϕ ja +1 ,b +1 − ϕ ja,b +1 (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | C ja,b = − M v ja,b +1 + ∆ t P ( (cid:28) ja,b ) (cid:16) ϕ ja,b − ϕ ja +1 ,b (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ϕ ja,b − ϕ ja +1 ,b +1 (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ϕ ja +1 ,b − ϕ ja +1 ,b +1 (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | D ja,b = − M v ja +1 ,b +1 + ∆ t P ( (cid:28) ja,b ) (cid:16) ϕ ja,b − ϕ ja +1 ,b (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ϕ ja,b +1 − ϕ ja,b (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ϕ ja,b +1 − ϕ ja +1 ,b (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | with n ( (cid:28) ja,b ) = ( ϕ ja +1 ,b − ϕ ja,b ) × ( ϕ ja,b +1 − ϕ ja,b ); n ( (cid:28) ja,b ) = ( ϕ a +1 ,b +1 − ϕ ja +1 ,b ) × ( ϕ ja,b − ϕ ja +1 ,b ) . n ( (cid:28) ja,b ) = ( ϕ a,b − ϕ ja,b +1 ) × ( ϕ ja +1 ,b +1 − ϕ ja,b +1 ); n ( (cid:28) ja,b ) = ( ϕ ja,b +1 − ϕ ja +1 ,b +1 ) × ( ϕ ja +1 ,b − ϕ ja +1 ,b +1 ) . mplicit integrator: The partial derivatives of the Lagrangian for a 2D barotropic fluidwith internal energy W ( ρ , J ), discretized under the mid-point rules, are given by A ja,b = ∆ t P ( (cid:28) ja,b ) (cid:16) ( ϕ ja +1 ,b + ϕ j +1 a +1 ,b ) − ( ϕ ja,b +1 + ϕ j +1 a,b +1 ) (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ( ϕ ja +1 ,b + ϕ j +1 a +1 ,b ) − ( ϕ ja +1 ,b +1 + ϕ j +1 a +1 ,b +1 ) (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ( ϕ ja +1 ,b +1 + ϕ j +1 a +1 ,b +1 ) − ( ϕ ja,b +1 + ϕ j +1 a,b +1 ) (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | B ja,b = ∆ t P ( (cid:28) ja,b ) (cid:16) ( ϕ ja,b +1 + ϕ j +1 a,b +1 ) − ( ϕ ja,b + ϕ j +1 a,b ) (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ( ϕ ja +1 ,b +1 + ϕ j +1 a +1 ,b +1 ) − ( ϕ ja,b + ϕ j +1 a,b ) (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ( ϕ ja +1 ,b +1 + ϕ j +1 a +1 ,b +1 ) − ( ϕ ja,b +1 + ϕ j +1 a,b +1 ) (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | C ja,b = ∆ t P ( (cid:28) ja,b ) (cid:16) ( ϕ ja,b + ϕ j +1 a,b ) − ( ϕ ja +1 ,b + ϕ j +1 a +1 ,b ) (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ( ϕ ja,b + ϕ j +1 a,b ) − ( ϕ ja +1 ,b +1 + ϕ j +1 a +1 ,b +1 ) (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ( ϕ ja +1 ,b + ϕ j +1 a +1 ,b ) − ( ϕ ja +1 ,b +1 + ϕ j +1 a +1 ,b +1 ) (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | D ja,b = ∆ t P ( (cid:28) ja,b ) (cid:16) ( ϕ ja,b + ϕ j +1 a,b ) − ( ϕ ja +1 ,b + ϕ j +1 a +1 ,b ) (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ( ϕ ja,b +1 + ϕ j +1 a,b +1 ) − ( ϕ ja,b + ϕ j +1 a,b ) (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + ∆ t P ( (cid:28) ja,b ) (cid:16) ( ϕ ja,b +1 + ϕ j +1 a,b +1 ) − ( ϕ ja +1 ,b + ϕ j +1 a +1 ,b ) (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | with n ( (cid:28) ja,b ) = (cid:16) ( ϕ ja +1 ,b + ϕ j +1 a +1 ,b ) − ( ϕ ja,b + ϕ j +1 a,b ) (cid:17) × (cid:16) ( ϕ ja,b +1 + ϕ j +1 a,b +1 ) − ( ϕ ja,b + ϕ j +1 a,b ) (cid:17) ; n ( (cid:28) ja,b ) = (cid:16) ( ϕ ja +1 ,b +1 + ϕ j +1 a +1 ,b +1 ) − ( ϕ ja +1 ,b + ϕ j +1 a +1 ,b ) (cid:17) × (cid:16) ( ϕ ja,b + ϕ j +1 a,b ) − ( ϕ ja +1 ,b + ϕ j +1 a +1 ,b ) (cid:17) . n ( (cid:28) ja,b ) = (cid:16) ( ϕ ja,b + ϕ j +1 a,b ) − ( ϕ ja,b +1 + ϕ j +1 a,b +1 ) (cid:17) × (cid:16) ( ϕ ja +1 ,b +1 + ϕ j +1 a +1 ,b +1 ) − ( ϕ ja,b +1 + ϕ j +1 a,b +1 ) (cid:17) ; n ( (cid:28) ja,b ) = (cid:16) ( ϕ ja,b +1 + ϕ j +1 a,b +1 ) − ( ϕ ja +1 ,b +1 + ϕ j +1 a +1 ,b +1 ) (cid:17) × (cid:16) ( ϕ ja +1 ,b + ϕ j +1 a +1 ,b ) − ( ϕ ja +1 ,b +1 + ϕ j +1 a +1 ,b +1 ) (cid:17) . Where the pressures P (cid:96) ( (cid:28) ja,b ), (cid:96) = 1 , , ,
4, were defined in (73).39 .2 Discrete Euler-Lagrange equations for 2D barotropic fluid ρ (cid:32) v ja,b − v j − a,b ∆ t (cid:33) = 14∆ s ∆ s (cid:40) − (cid:32) P ( (cid:28) ja,b ) (cid:16) ϕ ja,b +1 − ϕ ja +1 ,b (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + P ( (cid:28) ja − ,b ) (cid:16) ϕ ja − ,b − ϕ ja,b +1 (cid:17) × n ( (cid:28) ja − ,b ) | n ( (cid:28) ja − ,b ) | + P ( (cid:28) ja,b − ) (cid:16) ϕ ja +1 ,b − ϕ ja,b − (cid:17) × n ( (cid:28) ja,b − ) | n ( (cid:28) ja,b − ) | + P ( (cid:28) ja − ,b − ) (cid:16) ϕ ja,b − − ϕ ja − ,b (cid:17) × n ( (cid:28) ja − ,b − ) | n ( (cid:28) ja − ,b − ) | (cid:33) + (cid:32) P ( (cid:28) ja,b ) (cid:16) ϕ ja +1 ,b − ϕ ja +1 ,b +1 (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + P ( (cid:28) ja,b − ) (cid:16) ϕ ja +1 ,b − − ϕ ja +1 ,b (cid:17) × n ( (cid:28) ja,b − ) | n ( (cid:28) ja,b − ) | (cid:33) + (cid:32) P ( (cid:28) ja,b ) (cid:16) ϕ ja +1 ,b +1 − ϕ ja,b +1 (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | + P ( (cid:28) ja − ,b ) (cid:16) ϕ ja,b +1 − ϕ ja − ,b +1 (cid:17) × n ( (cid:28) ja − ,b ) | n ( (cid:28) ja − ,b ) | (cid:33) + (cid:32) P ( (cid:28) ja − ,b ) (cid:16) ϕ ja − ,b +1 − ϕ ja − ,b (cid:17) × n ( (cid:28) ja − ,b ) | n ( (cid:28) ja − ,b ) | + P ( (cid:28) ja − ,b − ) (cid:16) ϕ ja − ,b − ϕ ja − ,b − (cid:17) × n ( (cid:28) ja − ,b − ) | n ( (cid:28) ja − ,b − ) | (cid:33) + (cid:32) P ( (cid:28) ja − ,b − ) (cid:16) ϕ ja − ,b − − ϕ ja,b − (cid:17) × n ( (cid:28) ja − ,b − ) | n ( (cid:28) ja − ,b − ) | + P ( (cid:28) ja,b − ) (cid:16) ϕ ja,b − − ϕ ja +1 ,b − (cid:17) × n ( (cid:28) ja,b − ) | n ( (cid:28) ja,b − ) | (cid:33)(cid:41) Interpretation of the discrete Euler-Lagrange equations: discrete balance of momentum ρ ( v j − v j − ) / ∆ t = − [( P ext · L ext ) n ext − ( P int · L int ) n int ] j / area , where n ext , n int are unit vectors that point outward of the boundaries. φ a,bj φ d (a+1,b+1)(a,b) (a+1,b)(a+1,b-1)(a,b-1)(a-1,b-1)(a-1,b) (a,b+1)(a-1,b+1)1 23 4 1 23 41 23 41 23 4 φ a,b+1j φ a+1,b+1j φ a+1,bj φ a+1,b-1j φ a-1,b+1j φ a-1,bj φ a,b-1j φ a-1,b-1j L ext L int Figure 20:
Internal and external lengths L int , L ext . On the left: the force P ( (cid:28) ja,b ) (cid:16) ϕ ja +1 ,b − ϕ ja +1 ,b +1 (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | associated to the pressure P ( (cid:28) ja,b ) and the length | ϕ ja +1 ,b − ϕ ja +1 ,b +1 | . On the right P ( (cid:28) ja,b − ) (cid:16) ϕ ja +1 ,b − − ϕ ja +1 ,b (cid:17) × n ( (cid:28) ja,b − ) | n ( (cid:28) ja,b − ) | associated to the pressure P ( (cid:28) ja,b − ) and the length | ϕ ja +1 ,b − − ϕ ja +1 ,b | . Figure 22:
On the left: the force P ( (cid:28) ja,b ) (cid:16) ϕ ja,b +1 − ϕ ja +1 ,b (cid:17) × n ( (cid:28) ja,b ) | n ( (cid:28) ja,b ) | associated to the pressure P ( (cid:28) ja,b ) and the length | ϕ ja,b +1 − ϕ ja +1 ,b | . On the right P ( (cid:28) ja − ,b ) (cid:16) ϕ ja − ,b − ϕ ja,b +1 (cid:17) × n ( (cid:28) ja − ,b ) | n ( (cid:28) ja − ,b ) | associated to the pressure P ( (cid:28) ja − ,b ) and the length | ϕ ja − ,b − ϕ ja,b +1 | . A.3 Discrete Jacobian on (cid:28) ja,b,c in 3D J ( (cid:28) ja,b,c ) = ( F j a +1 ,b,c × F j a +1 ,b,c ) · F j a +1 ,b,c , J ( (cid:28) ja,b,c ) = ( F j a +1 ,b +1 ,c × F j a +1 ,b +1 ,c ) · F j a +1 ,b +1 ,c J ( (cid:28) ja,b,c ) = ( F j a,b +1 ,c × F j a,b +1 ,c ) · F j a,b +1 ,c , J ( (cid:28) ja,b,c ) = ( F j a,b +1 ,c +1 × F j a,b +1 ,c +1 ) · F j a,b +1 ,c +1 J ( (cid:28) ja,b,c ) = ( F j a,b,c +1 × F j a,b,c +1 ) · F j a,b,c +1 , J ( (cid:28) ja,b,c ) = ( F j a +1 ,b,c +1 × F j a +1 ,b,c +1 ) · F j a +1 ,b,c +1 J ( (cid:28) ja,b,c ) = ( F j a +1 ,b +1 ,c +1 × F j a +1 ,b +1 ,c +1 ) · F j a +1 ,b +1 ,c +1 . A.4 Derivatives of the discrete Lagrangian for 3D barotropic fluid
With the same definition as in § A.1, the partial derivatives of the discrete Lagrangianfor a 3D barotropic fluid with internal energy W ( ρ , J ) are listed below.41 ja,b,c = − M v ja,b,c + ∆ t (cid:110) − P (cid:16) n (12)1 + n (23)1 + n (31)1 (cid:17) + P n (12)4 + P n (23)3 + P n (31)2 (cid:111) ( (cid:28) ja,b,c ) ,B ja,b,c = − M v ja +1 ,b,c + ∆ t (cid:110) − P (cid:16) n (12)2 + n (23)2 + n (31)2 (cid:17) + P n (12)7 + P n (23)1 + P n (31)5 (cid:111) ( (cid:28) ja,b,c ) ,C ja,b,c = − M v ja,b +1 ,c + ∆ t (cid:110) − P (cid:16) n (12)3 + n (23)3 + n (31)3 (cid:17) + P n (12)6 + P n (23)5 + P n (31)1 (cid:111) ( (cid:28) ja,b,c ) ,D ja,b,c = − M v ja,b,c +1 + ∆ t (cid:110) − P (cid:16) n (12)4 + n (23)4 + n (31)4 (cid:17) + P n (12)1 + P n (23)7 + P n (31)6 (cid:111) ( (cid:28) ja,b,c ) ,E ja,b,c = − M v ja +1 ,b +1 ,c + ∆ t (cid:110) − P (cid:16) n (12)5 + n (23)5 + n (31)5 (cid:17) + P n (12)8 + P n (23)2 + P n (31)3 (cid:111) ( (cid:28) ja,b,c ) ,F ja,b,c = − M v ja,b +1 ,c +1 + ∆ t (cid:110) − P (cid:16) n (12)6 + n (23)6 + n (31)6 (cid:17) + P n (12)3 + P n (23)4 + P j n (31)8 (cid:111) ( (cid:28) ja,b,c ) ,G ja,b,c = − M v ja +1 ,b,c +1 + ∆ t (cid:110) − P (cid:16) n (12)7 + n (23)7 + n (31)7 (cid:17) + P n (12)2 + P n (23)8 + P n (31)4 (cid:111) ( (cid:28) ja,b,c ) ,H ja,b,c = − M v ja +1 ,b +1 ,c +1 + ∆ t (cid:110) − P (cid:16) n (12)8 + n (23)8 + n (31)8 (cid:17) + P n (12)5 + P n (23)6 + P n (31)7 (cid:111) ( (cid:28) ja,b,c ) , where n ( (cid:28) ja,b,c ) = (( ϕ ja +1 ,b,c − ϕ ja,b,c ) × ( ϕ ja,b +1 ,c − ϕ ja,b,c )) · ( ϕ ja,b,c +1 − ϕ ja,b,c ) , n ( (cid:28) ja,b,c ) = (( ϕ ja +1 ,b +1 ,c − ϕ ja +1 ,b,c ) × ( ϕ ja,b,c − ϕ ja +1 ,b,c )) · ( ϕ ja +1 ,b,c +1 − ϕ ja +1 ,b,c ) , n ( (cid:28) ja,b,c ) = (( ϕ ja,b,c − ϕ ja,b +1 ,c ) × ( ϕ ja +1 ,b +1 ,c − ϕ ja,b +1 ,c )) · ( ϕ ja,b +1 ,c +1 − ϕ ja,b +1 ,c ) , n ( (cid:28) ja,b,c ) = (( ϕ ja,b +1 ,c +1 − ϕ ja,b,c +1 ) × ( ϕ ja +1 ,b,c +1 − ϕ ja,b,c +1 )) · ( ϕ ja,b,c − ϕ ja,b,c +1 ) , n ( (cid:28) ja,b,c ) = (( ϕ ja,b +1 ,c − ϕ ja +1 ,b +1 ,c ) × ( ϕ ja +1 ,b,c − ϕ ja +1 ,b +1 ,c )) · ( ϕ ja +1 ,b +1 ,c +1 − ϕ ja +1 ,b +1 ,c ) , n ( (cid:28) ja,b,c ) = (( ϕ ja +1 ,b +1 ,c +1 − ϕ ja,b +1 ,c +1 ) × ( ϕ ja,b,c +1 − ϕ ja,b +1 ,c +1 )) · ( ϕ ja,b +1 ,c − ϕ ja,b +1 ,c +1 ) , n ( (cid:28) ja,b,c ) = (( ϕ ja,b,c +1 − ϕ ja +1 ,b,c +1 ) × ( ϕ ja +1 ,b +1 ,c +1 − ϕ ja +1 ,b,c +1 )) · ( ϕ ja +1 ,b,c − ϕ ja +1 ,b,c +1 ) , n ( (cid:28) ja,b,c ) = (( ϕ ja +1 ,b,c +1 − ϕ ja +1 ,b +1 ,c +1 ) × ( ϕ ja,b +1 ,c +1 − ϕ ja +1 ,b +1 ,c +1 )) · ( ϕ ja +1 ,b +1 ,c − ϕ ja +1 ,b +1 ,c +1 ) , and n (12)1 = ( ϕ ja +1 ,b,c − ϕ ja,b,c ) × ( ϕ ja,b +1 ,c − ϕ ja,b,c ) , n (23)1 = ( ϕ ja,b +1 ,c − ϕ ja,b,c ) × ( ϕ ja,b,c +1 − ϕ ja,b,c ) , n (31)1 = ( ϕ ja,b,c +1 − ϕ ja,b,c ) × ( ϕ ja +1 ,b,c − ϕ ja,b,c ) ,J ( (cid:28) ja,b,c ) = (( ϕ ja +1 ,b,c − ϕ ja,b,c ) × ( ϕ ja,b +1 ,c − ϕ ja,b,c )) · ( ϕ ja,b,c +1 − ϕ ja,b,c ) | s a +1 ,b,c − s a,b,c || s a,b +1 ,c − s a,b,c || s a,b,c +1 − s a,b,c | . Note that in n (12)1 , n (23)1 , n (31)1 we adopt the usual expressions (12), (23), (31) for permutation of theelements of the set { ( ϕ ja +1 ,b,c − ϕ ja,b,c ), ( ϕ ja,b +1 ,c − ϕ ja,b,c ), ( ϕ ja,b,c +1 − ϕ ja,b,c ) } which compose n . .5 Discrete Euler-Lagrange equations for 3D barotropic fluid ρ (cid:32) v ja,b,c − v j − a,b,c ∆ t (cid:33) = 18∆ s ∆ s ∆ s (cid:110) − P (cid:16) n (12)1 + n (23)1 + n (31)1 (cid:17) ( (cid:28) ja,b,c ) − P (cid:16) n (12)2 + n (23)2 + n (31)2 (cid:17) ( (cid:28) ja − ,b,c ) − P (cid:16) n (12)3 + n (23)3 + n (31)3 (cid:17) ( (cid:28) ja,b − ,c ) − P (cid:16) n (12)4 + n (23)4 + n (31)4 (cid:17) ( (cid:28) ja,b,c − ) − P (cid:16) n (12)5 + n (23)5 + n (31)5 (cid:17) ( (cid:28) ja − ,b − ,c ) − P (cid:16) n (12)6 + n (23)6 + n (31)6 (cid:17) ( (cid:28) ja,b − ,c − ) − P (cid:16) n (12)7 + n (23)7 + n (31)7 (cid:17) ( (cid:28) ja − ,b,c − ) − P (cid:16) n (12)8 + n (23)8 + n (31)8 (cid:17) ( (cid:28) ja − ,b − ,c − )+ (cid:16) P n (12)4 + P n (23)3 + P n (31)2 (cid:17) ( (cid:28) ja,b,c )+ (cid:16) P n (12)7 + P n (23)1 + P n (31)5 (cid:17) ( (cid:28) ja − ,b,c )+ (cid:16) P n (12)6 + P n (23)5 + P n (31)1 (cid:17) ( (cid:28) ja,b − ,c )+ (cid:16) P n (12)1 + P n (23)7 + P n (31)6 (cid:17) ( (cid:28) ja,b,c − )+ (cid:16) P n (12)8 + P n (23)2 + P n (31)3 (cid:17) ( (cid:28) ja − ,b − ,c )+ (cid:16) P n (12)3 + P n (23)4 + P n (31)8 (cid:17) ( (cid:28) ja,b − ,c − )+ (cid:16) P n (12)2 + P n (23)8 + P n (31)4 (cid:17) ( (cid:28) ja − ,b,c − )+ (cid:16) P n (12)5 + P n (23)6 + P n (31)7 (cid:17) ( (cid:28) ja − ,b − ,c − ) (cid:111) Interpretation of the discrete Euler-Lagrange equations: discrete balance of momentum ρ ( v j − v j − ) / ∆ t = − [( P ext · S ext ) n ext − ( P int · S int ) n int ] j / volume . (a+1,b-1)(a,b-1)(a-1,b-1) φ a+1,b-1j φ a,b-1j φ a-1,b-1j (a,b,c) (a+1,b,c)1 234 57 8 (a+1,b+1,c)(a+1,b+1,c+1)(a+1,b,c+1)(a,b+1,c+1)(a,b+1,c)(a,b,c+1) φ d φ a,b,cj φ a,b,c+1j φ a,b+1,c+1j φ a,b+1,cj φ a+1,b+1,c+1j φ a+1,b+1,cj φ a+1,b,c+1j φ a+1,b,cj F j F j F j F j F j F j Figure 23:
In red: the internal surfaces S int associated to the cell (cid:28) ja,b,c . References [1] Bauer, W. and Gay-Balmaz, F. [2019], Towards a geometric variational discretization of compress-ible fluids: the rotating shallow water equations,
J. Comp. Dyn. , (1), 1–37.[2] Bazaraa, M.S., Sherali, H.D., and Shetty, C.M. [2006] Nonlinear programming: Theory and Algo-rithms , Wiley, 2006.
3] Bridges, T. and Reich, S. [2001], Multi-symplectic integrators: numerical schemes for HamiltonianPDEs that conserve symplecticity,
Phys. Lett. A (4-5), 184–193.[4] Chorin, A.J. and Marsden, J. E. [1990],
A Mathematical Introduction to Fluid Dynamics . Springer.[5] Courant R. and Friedrichs K. O. [1948],
Supersonic Flow and Shock Waves , Institute for mathe-matics and mechanics, New York University, New-York, 1948.[6] Demoures, F., Gay-Balmaz, F., Desbrun, M., Ratiu, T. S., and Alejandro, A. [2017], A multisym-plectic integrator for elastodynamic frictionless impact problems,
Comput. Methods in Appl. Mech.Eng. , , 1025–1052.[7] Demoures, F., Gay-Balmaz, F., Kobilarov, M. and Ratiu T. S. [2014], Multisymplectic Lie groupvariational integrators for a geometrically exact beam in R , Commun. Nonlinear Sci. Numer.Simulat. , (10), 3492–3512.[8] Demoures, F., Gay-Balmaz, F., and Ratiu, T. S. [2014], Multisymplectic variational integrator andspace/time symplecticity, Anal. Appl. , (3), 341–391.[9] Demoures, F., Gay-Balmaz, F., and Ratiu, T. S. [2016], Multisymplectic variational integratorsfor nonsmooth Lagrangian continuum mechanics, Forum Math. Sigma, , e19, 54p.[10] Demoures, F. [2019], Multisymplectic variational integrators and constitutive discrete theory ofelasticity, submitted.[11] Donea J, Giuliani S, Halleux J. P. [1982], An arbitrary Lagrangian-Eulerian finite element methodfor transient dynamic fluid-structure interactions, Comput. Meth. in Appl. Mech. Eng. , , 689–723.[12] Donea, J. [1983], Arbitrary Lagrangian-Eulerian finite element methods. In T. Belytschko and T.J. R. Hughes, editors, Computational Methods in Transient Analysis , Elsevier.[13] Farhat, C., Rallu, A., Wang, K., and Belytschko, T. [2010], Robust and provably second-orderexplicit-explicit and implicit-explicit staggered time-integrators for highly non-linear compressiblefluid-structure interaction problems,
Int. J. Numer. Meth. Eng. , , 73–107.[14] Farkhutdinov, T., F. Gay-Balmaz, and V. Putkaradze [2020], Geometric variational approach tothe dynamics of porous media filled with incompressible fluid, Acta Mechanica , (9), 3897–3924. https://arxiv.org/pdf/2007.02605.pdf [15] Fetecau R.C., Marsden J.E., and West M. [2003], Variational multisymplectic formulations ofnonsmooth continuum mechanics, in Perspectives and Problems in Nonlinear Science , 229–261,Springer, New York, 2003.[16] Gawlik, E. S. and Gay-Balmaz, F. [2020], A variational finite element discretization of compressibleflow,
Found. Comput. Math. , 1-41. https://arxiv.org/pdf/1910.05648.pdf [17] Gay-Balmaz, F., Marsden, J. E. Ratiu, T. S. [2012], Reduced variational formulations in freeboundary continuum mechanics,
J. Nonlin. Sci. , (4), 463–497.[18] Gotay, M. J., J. Isenberg, J. E. Marsden, R. Montgomery, J. Sniatycki, P. B. Yasskin [1997],Momentum maps and classical fields. Part I: Covariant field theory, (1997), arXiv:physics/9801019v2 .[19] Gay-Balmaz, F. and V. Putkaradze [2020], Variational methods for fluid-structure interactions, in Springer Handbook of Variational Methods for Nonlinear Geometric Data , 175–205, Springer.[20] Hughes, T. J. R., Liu, W. K., Zimmerman, T. [1981], Lagrangian-Eulerian finite element formula-tion for incompressible viscous flow,
Comput. Meth. in Appl. Mech. Eng. , , 329–49.[21] Lew, A., Marsden, J. E., Ortiz, M., and West, M. [2003], Asynchronous variational integrators, Arch. Rational Mech. Anal. , (2), 85–146.[22] Lew, A., Marsden, J. E., Ortiz, M., and West, M. [2004], Variational time integrators, Internat. J.Numer. Methods Eng. , (1), 153–212.[23] Marsden, J. E. and Hughes, T. J. R. Mathematical Foundations of Elasticity . Prentice-Hall, 1983.
24] Marsden, J. E., Patrick, G. W., and Shkoller, S. [1998], Multisymplectic geometry, variationalintegrators and nonlinear PDEs,
Comm. Math. Phys. , , 351–395.[25] Marsden, J. E., Pekarsky, S. Shkoller, S., and West, M. [2001], Variational Methods, Multisym-plectic Geometry and Continuum Mechanics, J. Geom. Phys. , , 253-284.[26] Marsden, J. E. and West M. [2001], Discrete mechanics and variational integrators, Acta Numer. , 357–514.[27] Masud, A., Hughes, T. J. R. [1997], A space-time Galerkin/Least squares finite element formulationof the Navier-Stokes equations for moving domain problems, Comput. Meth. in Appl. Mech. Eng. , , 91–126.[28] Moreau, J.-J. [1973], On Unilateral Constraints, Friction and Plasticity , C.I.M.E. Summer Schools,1973.[29] Pavlov, D.
Structure-Preserving Discretization of Incompressible Fluids.
Thesis, Caltech, 2009.[30] Pavlov D., Mullen P., Tong Y., Kanso E., Marsden J. E., and Desbrun M. [2011], Structure-preserving discretization of incompressible fluids,
Physica D , (6), 443–458.[31] Rockafellar, R. T. [1970] Convex analysis , Princeton Univ. Press, Princeton.[32] Rockafellar, R. T. [1973]
Penalty methods and augmented Lagrangians in nonlinear programming ,in Fifth Conference on Optimization Techniques, R. Conti and A. Ruberti (eds.), Springer-Verlag,1973, 518–525.[33] Rockafellar, R. T. [1993], Lagrange multipliers and optimality,
SIAM Review (2), 183–238.[34] Rockafellar, R. T. and Wets, R. J.-B. [1998], Variational Analysis , Grundlehren der Mathematis-chen Wissenschaften, , Springer-Verlag, Berlin, 1998.[35] Soulaimani, A., Fortin, M., Dhatt, G., and Ouetlet, Y. [1991], Finite element simulation of two-and three-dimensional free surface flows,
Comput. Meth. in Appl. Mech. Eng. , , 265–296., 265–296.