Energy-Momentum Consistent Integration Scheme

• Energy-momentum consistent integration schemes are numerical methods that enforce the preservation or controlled decay of physical invariants such as energy and momentum at the discrete level.
• They utilize strategies like skew-symmetric formulations, discrete gradients, and variational principles to achieve long-term stability and accuracy in simulations.
• These methods are applied across disciplines including computational fluid dynamics, solid mechanics, and molecular dynamics to reliably simulate complex multiphysics systems.

An energy-momentum consistent integration scheme is a numerical methodology that guarantees the discrete conservation of (or a controlled, monotonic decay of) energy and momentum within a time-stepping framework, often even at the fully discrete level and despite the presence of nontrivial geometric, multiphysics, or high-Reynolds regimes. Such schemes are designed for dynamical partial differential equations (PDEs), ordinary differential equations (ODEs), or hybrid systems where the underlying physical model possesses conservation laws for energy and momentum (possibly including angular momentum). The class encompasses both implicit and explicit integrators, finite element and finite volume methods, as well as structure-preserving approaches for both continuum and particle dynamics problems. These methods are of fundamental importance in computational fluid dynamics, solid mechanics, molecular dynamics, phase-field modeling, and kinetic simulations where long-term fidelity and numerical stability of invariants or dissipative laws are required.

1. Principles of Energy-Momentum Consistent Integration

The key objective of energy-momentum consistent methods is the preservation (or correct dissipation) of physical invariants reflecting the first principles governing the system. The precise discrete formulation ensures that, at every integration step, quantities such as total kinetic energy, potential energy (or their combination), linear and angular momentum, and, when appropriate, additional invariants (Casimirs, enstrophy) obey exact or provably monotonic laws in complete analogy with their continuous counterparts.

These integrators utilize either skew-symmetric, split, or discrete-gradient discretizations of the relevant terms, and employ projection, Lagrange multiplier, or Galerkin constraints to strongly enforce conservation. In the presence of dissipation (viscosity, plasticity), discrete versions of the second law are realized by ensuring the entropy or energy decay rate is non-negative and matches the physical model.

In practice, conservation or dissipation at the discrete level is implemented via:

• Skew-symmetric or split advection terms to neutralize spurious numerical energy input,
• Upwind or central fluxes for mass, momentum, and energy with telescoping property in the discrete algebra,
• Carefully designed interfacial, penalty, or stabilization terms (e.g., interior penalty, flux jump penalizations) to avoid loss of conservation with equal-order interpolation,
• Discrete time-stepping algorithms (e.g., Crank–Nicolson, midpoint, discrete gradient) that are symmetric and allow exact telescoped conservation proofs,
• Variational, constrained, or directionality-based formulations (e.g., Livens’ principle, EMAC regularization, Dirac structure preservation).

2. Methodological Realizations Across Model Classes

The construction principles are concretely realized in a variety of settings, as follows:

Incompressible Flow and Hybrid Finite Element Methods

Hybrid finite element methods with discontinuous cell pressures and continuous facet-based flux multipliers enforce local and global mass and momentum balance as well as global discrete energy decay. The cellwise mass-flux and skew-symmetric upwinded advection structure combine with interior-penalty symmetric diffusive fluxes and a θ-scheme (θ ≥ ½ for energy stability), resulting in exact conservation of kinetic energy, local/global mass, and momentum at the discrete level—even in the absence of a solenoidal velocity field—without requiring a div-free reconstruction (Labeur et al., 2010).

Multiphase and Phase-Field Models

In two-phase flow, discrete energy-momentum consistency is achieved by including the phase-field correction to both the mass and momentum fluxes and using the same central-difference stencil for all transport terms and surface tension forces derived from the discrete free energy. This implementation allows for the exact preservation of mass, momentum, and (inviscid) kinetic energy, crucially permitting stable simulation at extreme density contrasts, and drastically reducing spurious interface-induced parasitic flows (Mirjalili et al., 2019).

Viscoelastodynamics and Thermomechanical Solids

Structure-preserving time-integrators, using midpoint-type balances and directionality constraints for the algorithmic stress (enforcing that discrete energy production by stress matches the discrete change in free energy and physical dissipation), yield fully discrete schemes that exactly enforce both energy balance and total (linear, angular) momentum conservation. The explicit update of internal inelastic variables satisfies this directionality, resulting in first- or second-order temporal schemes with exact dissipation matching and constraint satisfaction (Liu et al., 2023).

Hamiltonian Particle Dynamics and Molecular Systems

For discrete particle dynamics with pair, three-body, or many-body (EAM) potentials, energy-momentum integrators employ discrete gradient force assignments (Gonzalez, LaBudde–Greenspan), often requiring a coupled implicit solve per step, that provably conserve both total energy and linear momentum exactly. These methods are robust to mesh boundary conditions, support periodic systems, and avoid the drift and heating of standard explicit integrators (Schiebl et al., 2020, Marazzato et al., 2020). In large-scale $N$-body models, combinations of symmetric fast multipole solvers and hierarchically split time-steppers give exact pairwise force anti-symmetry and total momentum conservation at $O(N)$ complexity (Zhu, 2017).

Compressible Flows and Finite Volume Methods

Schemes that discretize the internal-energy equation (with corrective momentum–energy coupling terms) instead of the total-energy equation allow fully discrete conservation of the sum of kinetic and internal energy without resorting to Riemann solvers. The upwind flux construction, fractional-step pressure-correction, and appropriate source terms guarantee positivity and entropy bounds (Herbin et al., 2019).

Kinetic and Low-Rank Approximations

For kinetic systems such as the Vlasov–Poisson equation, dynamical low-rank approximations force the velocity-moment subspace to contain the exact moment basis, and the time-discretization is constructed to preserve the discrete continuity equations for mass, momentum, and energy exactly at every step (Einkemmer et al., 2021).

3. Variational and Structure-Preserving Discrete Mechanics

Advanced schemes utilize variational principles, notably the Livens (Hamilton–Pontryagin) action, to formulate mixed position–velocity–momentum updates in which the discrete conservation laws follow directly from the preservation properties of the discrete action. The use of Gonzalez discrete gradients ensures that the increment in energy is exactly matched by the increment in the discrete action, guaranteeing generalized energy preservation even in the presence of singular or semi-definite mass matrices, and extending to holonomically constrained and symmetric systems (Kinon et al., 2023).

4. Algorithmic Realizations and Stability Mechanisms

The specific schemes typically employ one or more of the following algorithmic strategies:

Scheme Type	Central Conservation Mechanism	Example Application
Skew-symmetric split (e.g., EMAC, upwind)	Nonlinear advection in skew-symmetric form	Incompressible flow (Labeur et al., 2010, Ingimarson, 2020)
Discrete gradient integrator (Labudde–Greenspan, Gonzalez)	Force derived from directionality (difference quotient or discrete gradient)	Hamiltonian and MD systems (Schiebl et al., 2020, Liu, 2022)
Explicit–synchronous or split time-stepping	Symmetric finite difference or explicit leapfrog with force integration	Hamiltonian ODEs (Marazzato et al., 2020), N-body (Zhu, 2017)
Midpoint or directionality constraint method	Midpoint (or higher-order) temporal discretization with stress/force update by matching free-energy increment	Viscoelasticity (Liu et al., 2023)
Finite volume with upwind, pressure correction	Positivity-preserving upwind convection together with total-energy correction	Compressible Euler (Herbin et al., 2019)

Energy stability and conservation require additional stabilization terms only for equal-order interpolation, non-divergence-constrained velocity fields, or in hyperbolic regimes. Choice of temporal weights (e.g., θ ≥ ½, χ=½) is necessary for energy decay; penalty parameters in finite element methods (e.g., interior penalty α) must be chosen “large enough” (typically $O(k^2)$, k = order of polynomial) to control associated interface jumps.

5. Analysis, Proofs, and Theoretical Guarantees

The core proofs rely on exploiting the telescoping property of the discrete variational forms, discrete summation-by-parts or flux-divergence theorems, and directionality (difference-quotient) identities for the energy. For most schemes, explicit test functions (e.g., $v=u_{n+\theta}$ for energy; $v=n_j$ for momentum) inserted into the discrete variational problem, combined with the skew-symmetry of the nonlinear terms and appropriateness of stabilization, yield exact discrete counterparts to the continuum conservation/dissipation laws. When the discrete structure matches the physical dissipation mechanism, the schemes deliver invariant preservation up to machine accuracy and strict monotonicity in the presence of physical damping.

Analysis results in the literature provide unconditional (modulo time-step for explicit variants) stability, proof of practical positivity (density, internal energy, pressure), and local or global discrete entropy inequalities, often with error estimates matching the theoretical order of the underlying spatial/temporal discretization (Labeur et al., 2010, Mirjalili et al., 2019, Liu et al., 2023, Herbin et al., 2019).

6. Implementation Strategy and Illustrative Applications

Efficient implementation leverages static condensation (eliminating cellwise variables in hybrid FE), projection or predictor-corrector splitting in velocity–pressure coupling, and parallel assembly of fluxes or force evaluations in molecular-dynamics applications. Automation via form compilers (e.g., UFL/FEniCS) permits the rapid translation of variational formulations into efficient low-level code (Labeur et al., 2010). Dynamical low-rank methods and tree-based N-body integrators exploit dimension reduction and fast summation for large-scale kinetic or gravitational systems.

Benchmark problems, including Kovasznay flow, backward-facing step (Navier–Stokes), chaotic advection, jet-in-crossflow (two-phase), mass-spring networks with singular $M$, and real-world molecular or planetary N-body simulations, robustly demonstrate the invariance properties, long-term stability, and accuracy of energy-momentum consistent integrators (Labeur et al., 2010, Mirjalili et al., 2019, Liu et al., 2023, Schiebl et al., 2020, Marazzato et al., 2020, Zhu, 2017, Liu, 2022, Einkemmer et al., 2021, Ingimarson, 2020).

7. Impact, Limitations, and Practical Recommendations

Energy-momentum consistent schemes have substantially advanced the simulation of multiphysics problems where fidelity of invariants is critical, enabling accurate coarse-mesh or large-step computations, supporting microcanonical accuracy in molecular ensembles, and ensuring stable long-time evolution free from artificial drift, energy blow-up, or mesh-induced degradation of momentum.

Limitations arise from the overhead of implicit solves (especially in discrete-gradient MD), potential loss of conservation with oversimplified stabilization, and challenges in preserving all invariants simultaneously in high-order or highly nonlinear PDEs without careful discretization design. Fine-tuning of penalty parameters and step sizes is sometimes essential, as is attention to machine-precision round-off effects for very long integrations.

Adoption of these integrators is currently standard in high-precision fluid dynamics, viscoelasticity, geophysical modeling, and molecular simulation, and is spreading in kinetic and multiphase flow applications as accuracy requirements and computational scale increase.

References

• "Energy stable and momentum conserving hybrid finite element method for the incompressible Navier-Stokes equations" (Labeur et al., 2010)
• "Consistent, energy-conserving momentum transport for simulations of two-phase flows using the phase field equations" (Mirjalili et al., 2019)
• "A continuum and computational framework for viscoelastodynamics: II. Strain-driven and energy-momentum consistent schemes" (Liu et al., 2023)
• "Energy-momentum conserving integration schemes for molecular dynamics" (Schiebl et al., 2020)
• "An explicit pseudo-energy conserving time-integration scheme for Hamiltonian dynamics" (Marazzato et al., 2020)
• "Consistent Internal Energy Based Schemes for the Compressible Euler Equations" (Herbin et al., 2019)
• "A consistent and conservative model and its scheme for $N$-phase-$M$-component incompressible flows" (Huang et al., 2021)
• "Energy-consistent integration of mechanical systems based on Livens principle" (Kinon et al., 2023)
• "A momentum conserving $N$-body scheme with individual timesteps" (Zhu, 2017)
• "On the design of energy-decaying momentum-conserving integrator for nonlinear dynamics using energy splitting and perturbation techniques" (Liu, 2022)
• "A mass, momentum, and energy conservative dynamical low-rank scheme for the Vlasov equation" (Einkemmer et al., 2021)
• "An energy, momentum and angular momentum conserving scheme for a regularization model of incompressible flow" (Ingimarson, 2020)