A canonical Hamiltonian formulation of the Navier–Stokes problem

Abstract

This paper presents a novel Hamiltonian formulation of the isotropic Navier–Stokes problem based on a minimum-action principle derived from the principle of least squares. This formulation uses the velocities $u_{i}(x_{j},t)$ and pressure $p(x_{j},t)$ as the field quantities to be varied, along with canonically conjugate momenta deduced from the analysis. From these, a conserved Hamiltonian functional $H^{*}$ satisfying Hamilton's canonical equations is constructed, and the associated Hamilton–Jacobi equation is formulated for both compressible and incompressible flows. This Hamilton–Jacobi equation reduces the problem of finding four separate field quantities ($u_{i}$,$p$) to that of finding a single scalar functional in those fields – Hamilton's principal functional ${S}^{*}[u_{i},p,t]$. Moreover, the transformation theory of Hamilton and Jacobi now provides a prescribed recipe for solving the Navier–Stokes problem: find ${S}^{*}$. If an analytical expression for ${S}^{*}$ can be obtained, it will lead via canonical transformation to a new set of fields which are simply equal to their initial values, giving analytical expressions for the original velocity and pressure fields. Failing that, if one can only show that a complete solution to this Hamilton–Jacobi equation does or does not exist, that will also resolve the question of existence of solutions. The method employed here is not specific to the Navier–Stokes problem or even to classical mechanics, and can be applied to any traditionally non-Hamiltonian problem.

1. Introduction

Given the title of this paper, it is incumbent on the authors to assure the reader that we do not claim to have done the impossible. A viscous fluid is, after all, a non-Hamiltonian system (Millikan 1929; Finlayson 1972a,b). There is no action integral for which Hamilton's principle (Hamilton 1833, 1834, 1835) yields the Navier–Stokes equations (Stokes 1845; Anderson, Tannehill & Pletcher 1984; Anderson 1995; Batchelor 2000; White 2006; Pozrikidis 2009; Cengel & Cimbala 2018) in their usual form (Millikan 1929; Finlayson 1972a,b), and we do not claim otherwise. Remarkably, however, a Hamiltonian formulation can still be found by considering a mathematically equivalent higher-order problem, as we will now demonstrate via simple example.

1.1. A motivating example

Consider the first-order initial-value problem

(1.1a,b)\begin{equation} \dot{v}={-}v, \quad v(0)=1, \end{equation}

with unique solution $v(t)={\rm e}^{-t}$. Here, $v(t)$ can be interpreted as the velocity of a lumped mass moving in a viscous medium in one dimension with linear damping. Like the traditional Navier–Stokes equations (Stokes 1845; Anderson et al. 1984; Anderson 1995; Batchelor 2000; White 2006; Pozrikidis 2009; Cengel & Cimbala 2018), this too is an intrinsically non-Hamiltonian problem, in that there is no action $\mathcal {S}$ for which Hamilton's principle ($\delta \mathcal {S}=0$) yields the governing equation $\dot {v}=-v$. Yet, if we simply differentiate both sides of the equation ($\ddot {v}=-\dot {v}$), use the original equation to write $\dot {v}=-v$ and apply the additional initial condition $\dot {v}(0)=-v(0)=-1$, we arrive at the mathematically equivalent second-order problem

(1.2a–c)\begin{equation} \ddot{v}=v, \quad v(0)=1, \quad \dot{v}(0)={-}1, \end{equation}

which has the same unique solution $v(t)={\rm e}^{-t}$ but which is Hamiltonian – not in the sense that the total mechanical energy is conserved, but in the sense that it has mathematically Hamiltonian structure.

As first observed by Sanders (2021, 2022, 2023a,b) and Sanders & Inman (2023), the associated action can be obtained by writing the original equation in standard form ($\mathcal {R}\equiv \dot {v}+v=0$), squaring the residual $\mathcal {R}$ and integrating over time

(1.3) \begin{align} \mathcal{S}^{*}[v]=\int\text{d}t\big(\tfrac{1}{2}\mathcal{R}^{2}\big)=\int \text{d}t\big[\tfrac{1}{2}\big(\dot{v}^{2}+2v\dot{v}+v^{2}\big)\big] \sim\int\text{d}t\big[\tfrac{1}{2}\big(\dot{v}^{2}+v^{2}\big)\big], \end{align}

where we have used the fact that $2v\dot {v}=\text {d}(v^{2})/\text {d}t$ is a total time derivative and can therefore be excluded from the action without changing the resulting Euler–Lagrange equation (Lanczos 1970). This is the so-called ‘time-averaged principle of least squares’ (Sanders 2021, 2022, 2023a,b; Sanders & Inman 2023): since $\mathcal {R}=0$ is a local minimum of $\mathcal {R}^{2}$, it is also a local minimum of $\int \text {d}t(\mathcal {R}^{2})$. Varying $v$, the first variation of $\mathcal {S}^{*}$ is

(1.4)\begin{equation} \delta\mathcal{S}^{*}=\int\text{d}t\left[\dot{v}\delta\dot{v}+v\delta v\right]=\int\text{d}t\left[(-\ddot{v}+v)\delta v\right]+\left[\dot{v}\delta v\right]_{t_{1}}^{t_{2}}, \end{equation}

yielding the second-order equation $\ddot {v}=v$ and revealing the canonically conjugate ‘momentum’ ${\rm \pi} \equiv \dot {v}$. Here, and in what follows, we will use the symbol ${\rm \pi}$ for canonically conjugate momenta, as is customary in Hamiltonian field theory, in order to avoid later confusion with the pressure field $p$. Since the mathematical constant $3.14159\ldots$ does not appear in the present work, there will be no ambiguity.

The corresponding Hamiltonian is obtained via the Legendre transform

(1.5) \begin{equation} H^{*}[v,{\rm \pi}]={\rm \pi}\dot{v}-\tfrac{1}{2}\big(\dot{v}^{2}+v^{2}\big)=\tfrac{1}{2}\big({\rm \pi}^{2}-v^{2}\big). \end{equation}

Notably, this Hamiltonian has nothing to do with the total mechanical energy of the system, although it is a conserved quantity. In fact, $H^{*}=0$ for the actual motion satisfying ${\rm \pi} \equiv \dot {v}=-v$. We note in passing that Liouville's theorem is satisfied, as the motion occurs along the line ${\rm \pi} =-v$, so that the phase-space volume, being always zero, is conserved. Hamilton's equations

(1.6a,b)\begin{equation} \dot{v}=\frac{\partial H^{*}}{\partial{\rm \pi}}, \quad \dot{\rm \pi}={-}\frac{\partial H^{*}}{\partial v}, \end{equation}

are mathematically equivalent to the second-order problem $\ddot {v}=v$ and therefore also mathematically equivalent to the original, first-order problem.

The associated Hamilton–Jacobi equation (Hamilton 1833, 1834, 1835; Jacobi 1837, 1842–1843; Whittaker 1904; Lanczos 1970) is

(1.7)\begin{equation} \frac{1}{2}\left(\frac{\partial{S}^{*}}{\partial v}\right)^{2}-\frac{1}{2}v^{2}+\frac{\partial{S}^{*}}{\partial t}=0, \end{equation}

where Hamilton's principal function ${S}^{*}={S}^{*}(v,t)$ serves as the generating function for a canonical transformation to a new coordinate $\phi$ which is constant and equal to its initial value. Although this is almost identical to the Hamilton–Jacobi equation for the simple harmonic oscillator – the only difference being the sign in front of $(1/2)v^{2}$ – the usual separable solution of the form ${S}^{*}(v,t)=W(v)+T(t)$, where $W(v)$ and $T(t)$ are functions of v and t (respectively), does not work, as the reader may check.

Instead, let us use a trial solution of the form

(1.8)\begin{equation} {S}^{*}(v,t)=F(t)v+\tfrac{1}{2}v^{2}+f(t), \end{equation}

where $F(t)$ and $f(t)$ are as yet undetermined functions of $t$. This trial solution was chosen to cancel the term $(1/2)v^{2}$ from the equation. Substituting our trial solution into the Hamilton–Jacobi equation, we find that

(1.9)\begin{equation} \tfrac{1}{2}[F(t)]^{2}+[F(t)+F'(t)]v+f'(t)=0. \end{equation}

In order for this equation to hold for all $v$, we must have the following:

(1.10)\begin{equation} F(t)+F'(t)=0 \quad \Rightarrow \quad F(t)=\alpha {\rm e}^{{-}t}, \end{equation}

where $\alpha$ is a constant of integration which will be used to transform to the new coordinate, and

(1.11)\begin{equation} \tfrac{1}{2}[F(t)]^{2}+f'(t)=0 \quad \Rightarrow \quad f(t)=\tfrac{1}{4}\alpha^{2}{\rm e}^{{-}2t}+\gamma, \end{equation}

where $\gamma$ is another constant of integration which is simply additive and can therefore be discarded.

In this way, we have that

(1.12)\begin{equation} {S}^{*}(v,t;\alpha)=\alpha {\rm e}^{{-}t}v+\tfrac{1}{2}v^{2}+\tfrac{1}{4}\alpha^{2}{\rm e}^{{-}2t}. \end{equation}

With one constant of integration ($\alpha$) to match the single degree of freedom ($v$), this is a complete solution to the Hamilton–Jacobi equation. The new coordinate $\phi$ (which is constant and equal to its initial value) is obtained via the canonical transformation

(1.13)\begin{equation} \phi=\frac{\partial{S}^{*}}{\partial\alpha}={\rm e}^{{-}t}v+\frac{1}{2}\alpha {\rm e}^{{-}2t}. \end{equation}

The numerical value of $\alpha$ is in turn obtained via the canonical relation

(1.14)\begin{equation} {\rm \pi}=\frac{\partial{S}^{*}}{\partial v}=\alpha {\rm e}^{{-}t}+v, \end{equation}

which, evaluated at $t=0$, gives $\alpha =-2$ (recall that ${\rm \pi} =\dot {v}$, and $\dot {v}(0)=-v(0)=-1$). Using the fact that the new coordinate $\phi$ is equal to its initial value, we have that

(1.15)\begin{equation} {\rm e}^{{-}t}v-{\rm e}^{{-}2t}=v(0)-1=0, \end{equation}

giving the correct solution $v(t)={\rm e}^{-t}$.

In summary, by doubling the order of the governing equation and supplying additional auxiliary conditions, we made a non-Hamiltonian problem into a Hamiltonian one (Sanders 2021, 2022, 2023a,b; Sanders & Inman 2023). Furthermore, this simple example demonstrates that the method correctly gives the solution to the original, non-Hamiltonian problem. Indeed, it would appear that every non-Hamiltonian problem belongs to an equivalence class of problems with the same solution, and within each such equivalence class there are Hamiltonian variants. The remainder of this paper applies that concept to the isotropic Navier–Stokes problem (Stokes 1845; Anderson et al. 1984; Anderson 1995; Batchelor 2000; White 2006; Pozrikidis 2009; Cengel & Cimbala 2018).

1.2. The Navier–Stokes problem

The incompressible Navier–Stokes equations (Stokes 1845; Anderson et al. 1984; Anderson 1995; Batchelor 2000; White 2006; Pozrikidis 2009; Cengel & Cimbala 2018) are given by

(1.16)$$\begin{gather} \rho \dot{u}_{i} + \rho u_{i,j}u_{j} + p_{,i} - \mu u_{i,jj} - \rho b_{i} = 0, \end{gather}$$

(1.17)$$\begin{gather}u_{i,i} = 0, \end{gather}$$

where $\rho$ is the constant and uniform density, $u_{i}=u_{i}(x_{j},t)$ is the velocity field, $p=p(x_{j},t)$ is the pressure field, $b_{i}=b_{i}(x_{j},t)$ is the body force field, subscript Roman indices label Euclidean tensor components ($i,j=1,2,3$), the $x_{j}$ are Eulerian spatial coordinates, $t$ is time, $\mu$ is the dynamic viscosity, a dot over a symbol denotes a partial time derivative ($\dot {u}_{i}=\partial {u}_{i}/\partial t$), a comma in a subscript indicates a spatial gradient ($p_{,i}=\partial p/\partial x_{i}$) and we employ the Einstein summation convention on repeated subscript indices. To be clear, the notation $u_{i}(x_{j},t)$ signifies that each component of the velocity field is a function of all three spatial coordinates $(x_{1},x_{2},x_{3})$ and time $t$. It has the same meaning as other common notations, such as $u_{i}(\boldsymbol {x},t)$ and $\boldsymbol {u}(\boldsymbol {x},t)$. Likewise for all other field quantities. In the case of a uniform gravitational field, $b_{i}=g_{i}$ coincides with the local acceleration due to gravity; however, in what follows, we make no assumptions about the functional form of $b_{i}(x_{j},t)$: it is completely arbitrary. There are four unknown field quantities: $u_{i}(x_{j},t)$ and $p(x_{j},t)$.

We seek, ultimately, a functional

(1.18)\begin{equation} H^{*}=H^{*}[u_{i},p,{\rm \pi}_{j},{\rm \pi}_{4};t], \end{equation}

where (${\rm \pi} _{i},{\rm \pi} _{4}$) are suitable ‘momenta’ conjugate to the field quantities ($u_{i},p$), such that Hamilton's canonical equations

(1.19a,b)$$\begin{gather} \dot{u}_{i}=\frac{\delta H^{*}}{\delta{\rm \pi}_{i}}, \quad\dot{p}=\frac{\delta H^{*}}{\delta{\rm \pi}_{4}}, \end{gather}$$

(1.20a,b)$$\begin{gather}\dot{\rm \pi}_{i}={-}\frac{\delta H^{*}}{\delta u_{i}}, \quad\dot{\rm \pi}_{4}={-}\frac{\delta H^{*}}{\delta p}, \end{gather}$$

constitute a mathematically equivalent second-order formulation of the problem, where $\delta H^{*}/\delta u_{i}$, $\delta H^{*}/\delta p$, $\delta H^{*}/\delta {\rm \pi}_{i}$ and $\delta H^{*}/\delta {\rm \pi}_{4}$ are the Volterra (1930) functional derivatives of $H^{*}$ with respect to the field quantities and the conjugate momenta. We will find that this is generally possible for a compressible fluid. For an incompressible fluid, the equation $\dot {p}=\delta H^{*}/\delta {\rm \pi}_{4}$ will need to be replaced by the incompressibility condition $u_{i,i}=0$, consistent with the well-known result that the pressure usually serves as a Lagrange multiplier for the incompressibility constraint (Lanczos 1970; Badin & Crisciani 2018) (refer to p. 361 of Lanczos (1970) and pp. 137 and 141 of Badin & Crisciani 2018).

The remainder of this paper is organized as follows. Section 2 gives a comprehensive overview of the relevant literature to date. Sections 3 and 4 contain the main results of the present work, culminating in a conserved Hamiltonian functional $H^{*}$ satisfying Hamilton's equations (1.19a,b) and (1.20a,b) for the mathematically equivalent second-order problem, along with the accompanying Hamilton–Jacobi equation (Hamilton 1833, 1834, 1835; Jacobi 1837, 1842–1843; Whittaker 1904; Lanczos 1970). Section 5 contains a discussion of the physical interpretation of the second-order formulation. Section 6 presents a brief case study in the form of one-dimensional flow over an infinite, flat plate. Finally, § 7 concludes the paper with a few closing remarks and an outline of how the present formulation can aid in resolving the question of existence and uniqueness of solutions to the Navier–Stokes problem.

By the end of the paper, we will have achieved precisely what the title promises: a canonical Hamiltonian formulation of the problem, opening new avenues toward resolution of one of the most famous unsolved problems in mathematics.

2. Literature review

The field of analytical mechanics, with foundations planted in Hamilton's principle of stationary action (Hamilton 1833, 1834, 1835) or d'Alembert's principle of virtual work (d'Alembert 1743), has been vital to the development of both classical and quantum physics since the eighteenth century. This approach is versatile and helpful to the physical understanding of the problem in question, and the foundation, structure and utility of Hamiltonian formalism is well documented (Becker 1954; Taylor 2005; Hamill 2014; Bohn 2018; Cline 2023; Fowler 2023). The supporting mathematics of the calculus of variations as well as symplectic and differential geometry can also be found in many excellent sources (Arnold 1989; Berndt 2001; Hall 2003; Boas 2006; Stone & Goldbart 2009; Gelfand & Fomin 2012; Arfken 2013; Needham 2021). It is therefore no surprise that researchers have been applying analytical formalism to classical fluids dating back to the time of Lagrange (Lagrange 1811; Lichtenstein 1929; Morrison 1998, 2006; Berdichevsky 2009a,b; dell'Isola & Gavrilyuk 2011; Badin & Crisciani 2018; Bedford 2021).

The task of obtaining solutions to the governing equations of fluid flow represents one of the most challenging problems in science and engineering. In most cases, the mathematical formulation is expressed as an initial-boundary-value problem: a set of coupled, nonlinear partial differential equations, which are to be solved subject to various initial and boundary conditions. The degree of complication of the governing equations depends on the type of the fluid. For a viscous fluid where the transport phenomena of friction and thermal conduction are included, the governing equations are called the Navier–Stokes equations (Stokes 1845; Anderson et al. 1984; Anderson 1995; Batchelor 2000; White 2006; Pozrikidis 2009; Cengel & Cimbala 2018). The Navier–Stokes equations are derived by applying fundamental physical principles – conservation of mass, conservation of momentum and conservation of energy – to a viscous fluid, and the derivation can be found in any fluid mechanics textbook (Anderson et al. 1984; Anderson 1995; Batchelor 2000; White 2006; Pozrikidis 2009; Cengel & Cimbala 2018). It is important to recognize that the Navier–Stokes equations as they are known today were not developed solely by Navier and Stokes; indeed, Poisson, Cauchy and others were also heavily involved in their development (Darrigol 2002). As far as the present authors are aware, to date there is still no firm answer to the question of whether or not there always exist unique, smooth, non-singular solutions to the three-dimensional Navier–Stokes equations (Lemarie-Rieusset 2018), and this constitutes one of the most famous unsolved problems in mathematics.

The application of analytical mechanics (Goldstein 1980; Arnold 1989; Fetter & Walecka 2003; Gelfand & Fomin 2012) to the field of fluid mechanics (Lanczos 1970) has recently seen a resurgence in interest (Salmon 1983, 1988; Brenier 2017; Giga, Kirshtein & Liu 2018; Mottaghi, Gabbai & Benaroya 2019; Taroco, Blanco & Feijoo 2020; Bedford 2021; Mavroeidis & Athanassoulis 2022) after a long history. In the absence of non-conservative forces, an inviscid fluid is a Hamiltonian system, and so the classical Hamiltonian theory applies. Serrin (1959), Benjamin (1984) and Holm, Marsden & Ratiu (1986) have all described variational and Hamiltonian formulations of incompressible, inviscid fluid flow. Roberts (1972) presented a Hamiltonian dynamic for weakly interacting vortices. This research obtained the canonical equations of Hamiltonian dynamics for a set of two well-separated vortex rings by setting up a Hamiltonian to define the set. Olver (1982) showed that the Euler equations of inviscid and incompressible fluid flow can be put into Hamiltonian form. Benjamin & Olver (1982) investigated the Hamiltonian structure of the water waves problem. They examined the symmetry groups of this problem, finding that Hamiltonian analysis enables the solution of conservative elements of the problem. However, the study also acknowledged that further study is needed to identify the physical significance of the mathematical results. Maddocks & Pego (1995) presented a novel Hamiltonian formulation of ideal fluid flow expressed in material coordinates. Their Hamiltonian formulation arises from a general approach for constrained systems that is not restricted to problems in fluid mechanics. Rather, it is widely applicable for obtaining unconstrained Hamiltonian dynamical systems from Lagrangian field equations that are subject to pointwise constraints. More recently, Arnold (2014) also studied the Hamiltonian nature of the ideal Euler equations.

Viscous forces are non-conservative, which presents a fundamental challenge when applying Hamilton's principle to viscous fluids (Millikan 1929; Finlayson 1972a,b; Lemarie-Rieusset 2018). Indeed, it is a well-known theorem (first proven by Millikan 1929) that the Navier–Stokes equations in their usual form cannot be derived from a classical action principle (Millikan 1929; Finlayson 1972a,b). Millikan (1929) summarizes his main result as follows:

It is impossible to derive the equations of steady motion of a viscous, incompressible fluid from a variation principle involving as Lagrangian function an expression in the velocity components and their first-order space derivatives, unless conditions are imposed on these velocity components such that all of the terms $vu_{,2}$, $wu_{,3}$, $wv_{,3}$, $uv_{,1}$, $uw_{,1}$, $vw_{,2}$ disappear from their positions in the Navier–Stokes equations (Millikan 1929).

(It should be noted that the six terms referred to above come from the convective acceleration $u_{i,j}u_{j}$, and Millikan (1929) uses the notation $u=u_{1}$, $v=u_{2}$ and $w=u_{3}$.) More generally, it has been shown that the existence of variational formulations is related to self-adjointness of the system with respect to a standard duality relation, a property that all non-conservative systems lack (Vainberg 1964). Within the last 80 years, many alternative methods have been developed in an attempt to circumvent the non-self-adjointness of dissipative systems (Prigogine & Glansdorff 1965; Biot 1970; Finlayson 1972a; Lebon & Lambermont 1973; Tonti 1973; Magri 1974; Telega 1979; Tonti 1984; Filippov 1989; Sieniutycz 2000; Robinson 2001; Galley 2013; Kim, Dargush & Lee 2016; Mottaghi et al. 2019; Taroco et al. 2020; Bersani & Caressa 2021; Junker & Balzani 2021). The mathematical study of alternative variational methods as applied to the Navier–Stokes equations in particular remains an ongoing endeavour (Oseledets 1989; Vujanovic & Jones 1989; Doering & Gibbon 1995; Fukagawa & Fujitani 2012; Jones 2015; Gay-Balmaz & Yoshimura 2017; Hieber, Robinson & Shibata 2017; Hochgerner 2018; Gay-Balmaz & Yoshimura 2019a,b; Rashad et al. 2021; Gonzalez & Taha 2022; Taha & Gonzalez 2022; Sanders 2023b).

Oseledets (1989) attempted to express the Navier–Stokes equations using Hamiltonian formalism. He was able to formalize the incompressible Euler equation but stated that his result is not valid for a compressible fluid. More recent attempts, such as Fukagawa & Fujitani (2012), Jones (2015) and Gay-Balmaz & Yoshimura (2017, 2019a,b), have enforced dissipation using a non-holonomic constraint on the entropy. Hochgerner (2018) attempted to obtain a Hamiltonian interacting particle system that could accurately model the fluid dynamics. His research separated the dynamics into slow (deterministic) and fast (stochastic) components to capture fine-scale effects. The study was able to derive the Navier–Stokes equation from a stochastic Hamiltonian system but ignored the stress tensor, was unable to separate configuration and momentum variables and did not establish energy conservation or dissipation.

Rashad et al. (2021) modelled the incompressible Navier–Stokes equations in so-called ‘port-Hamiltonian’ framework rather than the standard Hamiltonian framework. Their model used vector calculus instead of exterior calculus to minimize the number of operators. While the main goal of this research was increasing the interest of computational researchers in using vector calculus, they also demonstrated that vector calculus can help in the formulation of individual subsystems of Navier–Stokes equations and boundary ports of the model.

Gonzalez & Taha (2022), Taha & Gonzalez (2022) and Taha, Gonzalez & Shorbagy (2023) have recently applied Gauss's principle of least constraint (Gauss 1829) to the Navier–Stokes problem. Using Gauss's principle, Taha et al. (2023) have shown that, for an incompressible fluid, the magnitude of the pressure gradient is minimum over the domain, which they term the principle of minimum pressure gradient (PMPG). When applied to an inviscid fluid in two dimensions, the PMPG provides a closure condition based in first principles that could be considered a generalization of the Kutta condition to smooth geometries. It should be noted that Gauss's principle (Gauss 1829) is fundamentally different from Hamilton's principle (Hamilton 1833, 1834, 1835). Whereas the Hamiltonian framework involves an invariant action integral and employs variations in the coordinates (or, in continuum mechanics, the field quantities), Gauss's principle employs variations in the accelerations. As a result, the framework of Gauss's principle does not lead to canonical transformations.

Particularly relevant to the present work, Sanders (2021, 2022, 2023a,b) has shown that the higher-order dynamics are ‘intrinsically variational’, in the sense that higher-derivative versions of the classical equations of motion can be derived from a minimum-action principle even for dissipative systems, thus allowing inherently non-Hamiltonian problems to be treated as though they are Hamiltonian. This discovery has already led to two applications: the direct modal analysis of damped dynamical systems (Sanders 2022) and, subsequently, a new and more efficient algorithm for computing a damped system's resonant frequencies (Sanders & Inman 2023). Higher-derivative theories had been studied before in the realm of quantum gravity physics (Pais & Uhlenbeck 1950; Van den Berg & VanderVorst 2002; Kalies & VanderVorst 2004; Bender & Mannheim 2008; Smilga 2009; Mostafazadeh 2010; Baleanu et al. 2012; Chen et al. 2013) but until now they have not been applied to classical fluids. While the Navier–Stokes equations, in their standard form, may be unsuited to Hamiltonian formalism (Millikan 1929; Finlayson 1972a,b; Doering & Gibbon 1995; Hieber et al. 2017; Lemarie-Rieusset 2018), it will be shown here that the higher-order dynamics can be used to restate the problem in a form consistent with Hamiltonian and Hamilton–Jacobi formalism.

In conclusion, although the body of research surrounding the Navier–Stokes equations is extensive, it would appear that no canonical Hamiltonian formulation of the Navier–Stokes problem has been found to date. That is what the present work aims to achieve.

3. Lagrangian formulation of the problem

Although we are primarily interested in the incompressible form of the equations given by (1.16) and (1.17), here we will begin with the compressible form of the equations, with the understanding that we will eventually take the incompressible limit. For the compressible case, the linear momentum balance and continuity equations are given by

(3.1)$$\begin{gather} \mathcal{R}_{i}[u_{j},p,\rho;x_{k},t]\equiv \rho \dot{u}_{i} + \rho u_{i,j}u_{j} + p_{,i} - \mu u_{i,jj} - (\mu+\lambda)u_{j,ji} - \rho b_{i} = 0, \end{gather}$$

(3.2)$$\begin{gather}\mathcal{R}_{4}[u_{j},\rho]\equiv \dot{\rho} + \rho_{,i}u_{i} + \rho u_{i,i} = 0, \end{gather}$$

where $\rho =\rho (x_{j},t)$ is the density field (now one of the unknown field quantities along with $u_{i}$ and $p$), and $\lambda$ is an additional viscosity coefficient which, under Stokes's (Stokes 1845) hypothesis, is related to $\mu$ as $\lambda =-2\mu /3$, ensuring that the mechanical pressure agrees with the thermodynamic pressure. Henceforth, we will assume that all quantities have been suitably non-dimensionalized. The non-dimensional (constant) viscosities in (3.1) and (3.2) may be regarded as inverse Reynolds numbers, and the non-dimensional pressure may be considered to be normalized by the inertial scale $\rho _0 U^2$, with $\rho _0$ and $U$ appropriate density and velocity scales. As we will see later, starting from the compressible form of the equations will allow us to treat the pressure as a dynamical field variable alongside the velocities, rather than simply a Lagrange multiplier. Crucially, this will reveal in no uncertain terms what becomes of the momentum conjugate to the pressure (which will be identified later) in the incompressible limit.

In general, (3.1) and (3.2) would be appended with the energy equation, which introduces additional thermodynamic variables, such as temperature and enthalpy or entropy. Two of the thermodynamic variables are designated as ‘primary’, and equations of state are required to relate the remaining variables to the primary variables. Typically, pressure and temperature are chosen as the primary variables, and the equation of state for the density, for example, is expressed as $\rho =\rho (p,T)$. The conservation equations along with the equation of state constitute six equations for the six unknowns fields $(u_i,p,T,\rho )$. Henceforth in the present work, we will take the temperature to be constant, though we intend to consider variations in temperature in future work.

An incompressible flow is one for which the material derivative of the density vanishes, i.e. $\text {d}\rho /\text {d}t=\dot {\rho } + \rho _{,i}u_{i}=0$, and this condition serves as an equation of state. It is usually also assumed, for the sake of simplicity, that the density is both constant and uniform, further reducing the equation of state $\rho =\rho (p,T)$ to specification of $\rho =\rho _0$ as a system parameter. Consequently, (3.2) reduces to $\rho u_{i,i}=0$ and the energy equation is decoupled from the system. Accordingly, in the incompressible limit there are only four unknown field quantities $(u_{i},p)$ and the momentum balance and continuity equations suffice for the governing field equations.

We pause here to remark that all four field equations (3.1), (3.2) are first order in time with respect to the field quantities $u_{i}$ and $\rho$. This will be important shortly, when we double the order of the equations. It should also be noted, as mentioned previously, that the first-order problem described above is inherently non-Hamiltonian, in that there is no action $\mathcal {S}$ for which Hamilton's principle ($\delta \mathcal {S}=0$) yields the first-order field equations (Millikan 1929; Finlayson 1972a,b). Finally, we note that in the incompressible limit, $\mathcal {R}_{4}$ becomes independent of $\dot {\rho }$ and is no longer first order in time.

3.1. Second-order formulation

Although the first-order formulation of the problem is intrinsically non-Hamiltonian (Millikan 1929; Finlayson 1972a,b), nevertheless a Hamiltonian for the system may be found by considering a second-order formulation. Following Sanders (2023b), we observe that the actual motion of the fluid corresponds to the particular fields $(u_{i},p,\rho )$ for which the following action achieves a local minimum:

(3.3)\begin{equation} \mathcal{S}^{*}[u_{j},p,\rho]=\int\text{d}^{4}x(\tfrac{1}{2}\mathcal{R}_{i}\mathcal{R}_{i}+\tfrac{1}{2}\mathcal{R}_{4}\mathcal{R}_{4}), \end{equation}

where $\text {d}^{4}x=\text {d}\kern 0.06em x_{1}\,\text {d}\kern 0.06em x_{2}\,\text {d}\kern 0.06em x_{3}\,\text {d}t$, and the integral is carried out over both the control volume $\mathcal {V}$ occupied by the fluid ($x_{j}\in \mathcal {V}$) and the time period of interest ($t\in [t_{1},t_{2}]$). It must be emphasized that this action contains no new physics. Again, this is simply the principle of least squares (Finlayson 1972b) averaged over the space–time occupied by the fluid. The entire physics of the problem are already contained in the residuals ($\mathcal {R}_{i}$, $\mathcal {R}_{4}$).

Without an equation of state relating $\rho$ to $p$, the problem is underconstrained with five unknown field quantities and only four dynamical field equations. Anticipating the case of incompressible flow, where the density is constant and the four field quantities are $u_{i}$ and $p$, henceforth we will assume an equation of state of the form $\rho =\hat {\rho }(p)$, with $\hat {\rho }$ a known function determined either from first principles or empirically. In this way, the density field may be eliminated in favour of the pressure field, and the field equations assume the following form:

(3.4)$$\begin{gather} \mathcal{R}_{i}[u_{j},p;x_{k},t]\equiv \hat{\rho} \dot{u}_{i} + \hat{\rho} u_{i,j}u_{j} + p_{,i} - \mu u_{i,jj} - (\mu+\lambda)u_{j,ji} - \hat{\rho} b_{i} = 0, \end{gather}$$

(3.5)$$\begin{gather}\mathcal{R}_{4}[u_{j},p]\equiv \hat{\rho}'\dot{p} + \hat{\rho}'p_{,i}u_{i} + \hat{\rho} u_{i,i} = 0, \end{gather}$$

where $\hat {\rho }'=\text {d}\hat {\rho }/\text {d}p$. We note that, under equilibrium conditions, $\hat {\rho }'$ is related to the speed of sound $c$ and the bulk modulus $K$ as $\hat {\rho }'=1/c^{2}=\rho /K$ (for incompressible fluids, $\hat {\rho }'\equiv 0$ and the speed of sound and bulk modulus are both infinite). Having specified $\hat {\rho }(p)$, $\mu$, $\lambda$ and $b_{i}(x_{j},t)$, and having prescribed appropriate auxiliary conditions (initial and boundary conditions), one seeks the four field quantities $(u_{i},p)$ satisfying the governing field equations and the auxiliary conditions. To recover the case of incompressible flow, we will eventually take $\hat {\rho }'\equiv 0$.

We pause here to note that, even though the residuals ($\mathcal {R}_{i}$, $\mathcal {R}_{4}$) vanish for the actual motion, they are not trivially zero. That is, the residuals only vanish for the particular fields $(u_{i},p)$ which satisfy the first-order field equations (3.4) and (3.5); they do not vanish for every conceivable $(u_{i},p)$. Thus it is not appropriate to take $\mathcal {R}_{i}\equiv 0$, $\mathcal {R}_{4}\equiv 0$. We will return to this point later when we discuss the Hamiltonian formulation of the problem.

For now, we note that the action $\mathcal {S}^{*}=\mathcal {S}^{*}[u_{i},p]$ defines a Lagrangian

(3.6)\begin{equation} L^{*}[u_{i},p;t]=\int \text{d}^{3}x(\mathcal{L}^{*}) , \end{equation}

where the integral is carried out over the volume $\mathcal {V}$ only ($\text {d}^{3}x=\text {d}\kern 0.06em x_{1}\,\text {d}\kern 0.06em x_{2}\,\text {d}\kern 0.06em x_{3}$), with Lagrangian density

(3.7)\begin{equation} \mathcal{L}^{*}[u_{j},p;x_{k},t]=\tfrac{1}{2}\mathcal{R}_{i}\mathcal{R}_{i}+\tfrac{1}{2}\mathcal{R}_{4}\mathcal{R}_{4}. \end{equation}

Because the residuals ($\mathcal {R}_{i}, \mathcal {R}_{4}$) have been non-dimensionalized, the Lagrangian density is also dimensionless. Once again, even though the Lagrangian vanishes for the actual motion, it is not trivially zero, and it is not appropriate to take $L^{*}\equiv 0$.

As noted above, the actual motion of the fluid corresponds to the particular fields $(u_{i},p)$ for which $\mathcal {S}^{*}$ achieves a local minimum. To obtain the Euler–Lagrange equations, the conjugate momenta and the natural auxiliary conditions, we insist that $\mathcal {S}^{*}$ not vary to first order ($\delta \mathcal {S}^{*}=0$) under small variations in the fields $(\delta u_{i},\delta p)$. Evaluating $\delta \mathcal {S}^{*}$, integrating by parts, and collecting like terms, we find that

(3.8a)\begin{align} \delta\mathcal{S}^{*}&=\int\text{d}^{4}x\left\{\left[-\frac{\partial}{\partial t} \left(\hat{\rho}\mathcal{R}_{i}\right) -\frac{\partial}{\partial x_{j}} \left(\hat{\rho} \mathcal{R}_{i}u_{j}\right)+ \hat{\rho} \mathcal{R}_{j}u_{j,i}- \mu \mathcal{R}_{i,jj} \right.\right. \end{align}

(3.8b)\begin{align} &\quad - \left.(\mu+\lambda)\mathcal{R}_{j,ij}+ \hat{\rho}'\mathcal{R}_{4}p_{,i}-\frac{\partial}{\partial x_{i}}\left( \hat{\rho} \mathcal{R}_{4}\right)\right]\delta{u}_{i} \end{align}

(3.8c)\begin{align} &\quad +\left[\vphantom{\frac{\partial}{\partial t}}\hat{\rho}'\mathcal{R}_{i} \dot{u}_{i}+ \hat{\rho}' \mathcal{R}_{i}u_{i,j}u_{j}- \mathcal{R}_{i,i}- \hat{\rho}' \mathcal{R}_{i}b_{i} + \hat{\rho}''\mathcal{R}_{4}\dot{p}\right. \end{align}

(3.8d)\begin{align} &\quad -\left.\left.\frac{\partial}{\partial t}\left( \hat{\rho}'\mathcal{R}_{4}\right) + \hat{\rho}''\mathcal{R}_{4}p_{,i}u_{i} - \frac{\partial}{\partial x_{i}}\left( \hat{\rho}'\mathcal{R}_{4}u_{i}\right)+ \hat{\rho}' \mathcal{R}_{4}u_{i,i}\right] \delta p \right\} \end{align}

(3.8e)\begin{align} &\quad +\int\text{d}^{3}x\left[\hat{\rho}\mathcal{R}_{i} \delta{u}_{i}+\hat{\rho}'\mathcal{R}_{4}\delta{p}\right]_{t_{1}}^ {t_{2}} \end{align}

(3.8f)\begin{align} &\quad +\int\text{d}^{2}x\,\text{d}t\left\{\left[\hat{\rho} \mathcal{R}_{i}u_{j}n_{j}+\mu \mathcal{R}_{i,j}n_{j}+ (\mu+\lambda)\mathcal{R}_{j,i}n_{j}+\hat{\rho} \mathcal{R}_{4}n_{i}\right]\delta u_{i}\right. \end{align}

(3.8g)\begin{align} &\quad +\left.\left[- \mu \mathcal{R}_{i}n_{j}- (\mu+\lambda)\mathcal{R}_{j}n_{i}\right]\delta u_{i,j}+\left[\mathcal{R}_{i}n_{i}+\hat{\rho}'\mathcal{R}_{4}u_{i}n_{i}\right]\delta p\right\}, \end{align}

where the purely volumetric integral ($\text {d}^{3}x$) is carried out over $\mathcal {V}$, the surface integral ($\text {d}^{2}x$) is carried out over the boundary $\partial \mathcal {V}$ and $n_{i}$ is the unit outward normal vector to $\partial \mathcal {V}$. Note that, because we are using Eulerian coordinates $x_{j}$, the volume element $\text {d}^{3}x$ is not to be varied.

The Euler–Lagrange equations (which hold for all $x_{j}\in \mathcal {V}$) may be read directly from the space–time ($\text {d}^{4}x$) integral

(3.9)\begin{align} \delta{u}_{i}: \quad &{-}\frac{\partial}{\partial t} \left(\hat{\rho}\mathcal{R}_{i}\right) -\frac{\partial}{\partial x_{j}} \left(\hat{\rho} \mathcal{R}_{i}u_{j}\right)+ \hat{\rho} \mathcal{R}_{j}u_{j,i}- \mu \mathcal{R}_{i,jj} \nonumber\\ &- (\mu+\lambda)\mathcal{R}_{j,ij}+ \hat{\rho}'\mathcal{R}_{4}p_{,i}-\frac{\partial}{\partial x_{i}}\left( \hat{\rho} \mathcal{R}_{4}\right)=0, \end{align}

(3.10)\begin{align} \delta p:\quad & \hat{\rho}'\mathcal{R}_{i} \dot{u}_{i}+ \hat{\rho}' \mathcal{R}_{i}u_{i,j}u_{j}- \mathcal{R}_{i,i}- \hat{\rho}' \mathcal{R}_{i}b_{i} + \hat{\rho}''\mathcal{R}_{4}\dot{p} \nonumber\\ &-\frac{\partial}{\partial t}\left( \hat{\rho}'\mathcal{R}_{4}\right) + \hat{\rho}''\mathcal{R}_{4}p_{,i}u_{i} - \frac{\partial}{\partial x_{i}}\left(\hat{\rho}'\mathcal{R}_{4}u_{i}\right)+ \hat{\rho}' \mathcal{R}_{4}u_{i,i} =0. \end{align}

It should be noted that all four Euler–Lagrange equations (3.9), (3.10) are second order in time, as they involve time derivatives of the residuals. By doubling the order of the equations, we have put the problem in Hamiltonian form, consistent with the general result of Sanders (2023b). We also note that all four Euler–Lagrange equations of the second-order formulation are automatically satisfied by the solution to the first-order formulation (i.e. the actual motion), for which $\mathcal {R}_{i}=0$ and $\mathcal {R}_{4}=0$ everywhere and at all times.

Corresponding to each of the four field quantities is a canonically conjugate ‘momentum’ field, which can be read from (3.8e). The momenta conjugate to the velocities $u_{i}$ are

(3.11)\begin{equation} {\rm \pi}_{i}\equiv \hat{\rho}\mathcal{R}_{i}, \end{equation}

and the momentum conjugate to the pressure $p$ is

(3.12)\begin{equation} {\rm \pi}_{4}\equiv \hat{\rho}'\mathcal{R}_{4}. \end{equation}

In the forthcoming Hamiltonian formulation, the conjugate momenta will be used to eliminate the (partial) time derivatives $(\dot {u}_{i},\dot {p})$ of the field quantities from the Hamiltonian. In general, Hamilton's principle would insist that the variations $(\delta u_{i},\delta p)$ vanish at the endpoints $t=t_{1}$ and $t=t_{2}$ to ensure that the purely volumetric integral (3.8e) vanishes identically. Interestingly, for the actual motion ($\mathcal {R}_{i}=0,\mathcal {R}_{4}=0$), the volumetric integral (3.8e) already vanishes even without taking $(\delta u_{i},\delta p)$ to vanish at $t_{1}$ and $t_{2}$. We interpret this to mean that the actual motion is the natural evolution of the second-order formulation (Sanders 2023b).

Although the conjugate momenta (${\rm \pi} _{i}$, ${\rm \pi} _{4}$) do not coincide with conventional linear or angular momenta, there is nonetheless a curious mathematical connection between the conjugate momenta and the linear momentum density $P_{i}=\rho u_{i}$, which we will see presently from the natural boundary conditions. These are read directly from the surface ($\text {d}^{2}x$) integral

(3.13)$$\begin{gather} \delta u_{i}:\quad \hat{\rho} \mathcal{R}_{i}u_{j}n_{j}+\mu \mathcal{R}_{i,j}n_{j}+ (\mu+\lambda)\mathcal{R}_{j,i}n_{j}+\hat{\rho} \mathcal{R}_{4}n_{i}=0, \end{gather}$$

(3.14)$$\begin{gather}\delta u_{i,j}:\quad - \mu \mathcal{R}_{i}n_{j}- (\mu+\lambda)\mathcal{R}_{j}n_{i}=0, \end{gather}$$

(3.15)$$\begin{gather}\delta p:\quad\mathcal{R}_{i}n_{i}+\hat{\rho}'\mathcal{R}_{4}u_{i}n_{i}=0. \end{gather}$$

This last condition, (3.15), establishes a connection between the new conjugate momenta and the conventional linear momenta. Multiplying (3.15) by $\hat {\rho }$, and noting that $\hat {\rho }u_{i}=\rho u_{i}= P_{i}$, we find that

(3.16)\begin{equation} ({\rm \pi}_{i}+{\rm \pi}_{4}P_{i})n_{i}=0. \end{equation}

Evidently, boundary condition (3.15) states that the flux of the vector $\varPi _{i}\equiv {\rm \pi}_{i}+{\rm \pi} _{4}P_{i}$ through the boundary $\partial \mathcal {V}$ should vanish. It is interesting that this new vector $\varPi _{i}$ contains both old and new momenta, with ${\rm \pi} _{4}$ ‘carried’ (i.e. given direction) by $P_{i}$. The actual physical meaning of these natural boundary conditions is less clear and may require further investigation.

3.2. Equivalence of the first- and second-order formulations

The first- and second-order formulations are mathematically equivalent, in the sense that imposing identical auxiliary conditions on the two formulations will yield identical solutions $(u_{i},p)$. In other words, with identical auxiliary conditions, $(u_{i},p)$ is a solution to the first-order formulation if and only if the same $(u_{i},p)$ is a solution to the second-order formulation.

The proof is straightforward. Consider the two formulations separately, and impose on the second-order formulation identical auxiliary conditions to those of the first-order formulation. In particular, just like the simple example given in § 1.1, the second-order formulation requires additional auxiliary conditions over and above those applied to the first-order formulation. These include initial conditions making $\mathcal {R}_{i}(x_{k},0)=0$ and $\mathcal {R}_{4}(x_{k},0)=0$ for all $x_{k}\in \mathcal {V}\cup \partial \mathcal {V}$, along with boundary conditions making $\mathcal {R}_{i}(x_{k},t)=0$, $\mathcal {R}_{i,j}(x_{k},t)=0$ and $\mathcal {R}_{4}(x_{k},t)=0$ for all $x_{k}\in \partial \mathcal {V}$ and all times $t$. By supposition, the auxiliary conditions applied to the two formulations are identical, so it suffices to show that $(u_{i},p)$ satisfies the governing field equations (3.4), (3.5) of the first-order formulation ($\mathcal {R}_{i}=0$ and $\mathcal {R}_{4}=0$) everywhere in $\mathcal {V}$ and at all times if and only if $(u_{i},p)$ satisfies the Euler–Lagrange equations (3.9), (3.10) of the second-order formulation everywhere in $\mathcal {V}$ and at all times.

Suppose first that $(u_{i},p)$ satisfies the governing field equations (3.4), (3.5) of the first-order formulation everywhere in $\mathcal {V}$ and at all times. Then $\mathcal {R}_{i}=0$, $\mathcal {R}_{4}=0$, and $(u_{i},p)$ is a trivial solution to the Euler–Lagrange equations (3.9), (3.10) of the second-order formulation.

Conversely, suppose that $(u_{i},p)$ satisfies the Euler–Lagrange equations (3.9), (3.10) of the second-order formulation everywhere in $\mathcal {V}$ and at all times. We note that $(\mathcal {R}_{i}=0,\mathcal {R}_{4}=0)$ constitutes an equilibrium solution of the Euler–Lagrange equations (3.9), (3.10). Thus, because the initial conditions are chosen such that $\mathcal {R}_{i}(x_{k},0)=0$ and $\mathcal {R}_{4}(x_{k},0)=0$ for all $x_{k}\in \mathcal {V}\cup \partial \mathcal {V}$, and because the boundary conditions are chosen such that $\mathcal {R}_{i}(x_{k},t)=0$, $\mathcal {R}_{i,j}(x_{k},t)=0$ and $\mathcal {R}_{4}(x_{k},t)=0$ for all $x_{k}\in \partial \mathcal {V}$ and all times $t$, then $\mathcal {R}_{i}$ and $\mathcal {R}_{4}$ will remain identically zero everywhere in $\mathcal {V}$ for all future times. Thus, $(u_{i},p)$ satisfies the governing field equations (3.4), (3.5) of the first-order formulation everywhere in $\mathcal {V}$ and at all times. This completes the proof, and we have established that the two formulations are equivalent.

4. Hamiltonian formulation of the problem

We are now ready to proceed with the Hamiltonian formulation of the problem. The Lagrangian $L^{*}$ has a corresponding Hamiltonian

(4.1)\begin{equation} H^{*}=\int \text{d}^{3}x(\mathcal{H}^{*}) , \end{equation}

with the Hamiltonian density $\mathcal {H}^{*}$ obtained from the Lagrangian density $\mathcal {L}^{*}$ via the Legendre transform

(4.2)\begin{equation} \mathcal{H}^{*}={\rm \pi}_{i}\dot{u}_{i}+{\rm \pi}_{4}\dot{p}-\mathcal{L}^{*}={\rm \pi}_{i}\dot{u}_{i}+{\rm \pi}_{4}\dot{p}-\tfrac{1}{2}\mathcal{R}_{i}\mathcal{R}_{i}-\tfrac{1}{2}\mathcal{R}_{4}\mathcal{R}_{4}. \end{equation}

Again, this $H^{*}$ has nothing to do with the total mechanical energy of the system, although it is a conserved quantity, since $H^{*}=0$ for the actual motion – just as in the example of § 1.1. In order to write down Hamilton's equations, we must express $\mathcal {H}^{*}$ in terms of the fields and the conjugate momenta, eliminating $\dot {u}_{i}$ and $\dot {p}$.

Observe that $\mathcal {R}_{i}={\rm \pi} _{i}/\hat {\rho }$, and ignoring for the moment the incompressible limit, we may write $\mathcal {R}_{4}={\rm \pi} _{4}/\hat {\rho }'$ ($\hat {\rho }'\neq 0$). In this way, using the functional expressions for the residuals given by (3.4) and (3.5), we find that

(4.3)$$\begin{gather} \dot{u}_{i}=\frac{1}{(\hat{\rho})^{2}}{{\rm \pi}_{i}}-\frac{1}{\hat{\rho}}\left(\hat{\rho}u_{i,j}u_{j} +{ p_{,i} - \mu u_{i,jj} - (\mu+\lambda)u_{j,ji}} - \hat{\rho}b_{i}\right); \end{gather}$$

(4.4)$$\begin{gather}\dot{p}=\frac{1}{(\hat{\rho}')^{2}}{{\rm \pi}_{4}}-\frac{1}{\hat{\rho}'}\left({\hat{\rho}'p_{,i}u_{i} + \hat{\rho}u_{i,i}}\right), \quad \hat{\rho}'\neq0; \end{gather}$$

and

(4.5)\begin{align} \mathcal{H}^{*}[u_{i},p,{\rm \pi}_{j},{\rm \pi}_{4};x_{k},t]&=\frac{1}{2}\frac{1}{(\hat{\rho})^{2}}{{\rm \pi}_{i}{\rm \pi}_{i}}-\frac{1}{\hat{\rho}}(\hat{\rho}u_{i,j}u_{j} +{p_{,i} - \mu u_{i,jj} - (\mu+\lambda)u_{j,ji}} - \hat{\rho}b_{i}){\rm \pi}_{i}\nonumber\\ &\quad +\frac{1}{2}\frac{1}{(\hat{\rho}')^{2}}{{\rm \pi}_{4}{\rm \pi}_{4}}-\frac{1}{\hat{\rho}'}({\hat{\rho}'p_{,i}u_{i} + \hat{\rho}u_{i,i}}){\rm \pi}_{4}, \quad \hat{\rho}'\neq0. \end{align}

Hamilton's equations (Hamilton 1834, 1835) are as follows:

(4.6a,b)$$\begin{gather} \dot{u}_{i}=\frac{\delta H^{*}}{\delta{\rm \pi}_{i}}, \quad\dot{p}=\frac{\delta H^{*}}{\delta{\rm \pi}_{4}}, \end{gather}$$

(4.7a,b)$$\begin{gather}\dot{\rm \pi}_{i}={-}\frac{\delta H^{*}}{\delta u_{i}}, \quad \dot{\rm \pi}_{4}={-}\frac{\delta H^{*}}{\delta p}, \end{gather}$$

where $\delta H^{*}/\delta u_{i}$, $\delta H^{*}/\delta p$, $\delta H^{*}/\delta {\rm \pi}_{i}$ and $\delta H^{*}/\delta {\rm \pi}_{4}$ are the Volterra (Volterra 1930) functional derivatives of $H^{*}$ with respect to the field quantities and the conjugate momenta. Equations (4.6a,b) reproduce (4.3) and (4.4), respectively. Equations (4.7a,b) in turn reproduce the Euler–Lagrange equations (3.9), (3.10) of the second-order formulation.

We return now to our previous observation concerning the vanishing of the residuals. While $H^{*}$ vanishes for the particular fields $(u_{i},p)$ that satisfy the governing field equations (3.4), (3.5) of the first-order formulation, it does not vanish for every conceivable $(u_{i},p)$. The latter would imply, according to (4.6a,b), that $\dot {u}_{i}\equiv 0$ and $\dot {p}\equiv 0$, which is not generally the case. This observation, and the fact that (4.7a,b) faithfully reproduce the Euler–Lagrange equations (3.9), (3.10) of the second-order formulation, confirm that the Hamiltonian formulation described above is, in fact, a legitimate reformulation of the problem. In the following section, we develop the Hamilton–Jacobi theory as it relates to the present formulation, the goal being to find a canonical transformation to a new set of fields ($\phi _{i}, \phi _{4}$) and conjugate momenta for which the Hamiltonian does vanish identically.

To obtain the Hamiltonian for incompressible flow, we set $\hat {\rho }'\equiv 0$ from the beginning, in which case $\mathcal {R}_{4}$ reduces to $\hat {\rho } u_{i,i}$ and ${\rm \pi} _{4}$ vanishes identically, consistent with the fact that $\mathcal {R}_{4}$ becomes independent of $\dot {p}$. The Hamiltonian density in turn reduces to

(4.8)\begin{equation} \mathcal{H}^{*}={\rm \pi}_{i}\dot{u}_{i}-\tfrac{1}{2}\mathcal{R}_{i}\mathcal{R}_{i}, \quad \hat{\rho}'\equiv0, \end{equation}

or, in terms of the conjugate momenta,

(4.9)\begin{equation} \mathcal{H}^{*}[u_{i},p,{\rm \pi}_{j};x_{k},t]=\frac{1}{2}\frac{1}{{\rho}^{2}}{{\rm \pi}_{i}{\rm \pi}_{i}}-\frac{1}{\rho}({\rho}u_{i,j}u_{j} + p_{,i} - \mu u_{i,jj} - {\rho}b_{i}){\rm \pi}_{i}, \quad \hat{\rho}'\equiv0,\end{equation}

where the density $\hat {\rho }=\rho$ is a constant and we have used the fact that $u_{i,i}=0$. Hamilton's equations $\dot {u}_{i}=\delta H^{*}/\delta {\rm \pi}_{i}$, $\dot {{\rm \pi} }_{i}=-\delta H^{*}/\delta u_{i}$ and $0\equiv \dot {{\rm \pi} }_{4}=-\delta H^{*}/\delta p$ still apply (and reproduce the corresponding equations in the incompressible limit), but $\dot {p}=\delta H^{*}/\delta {\rm \pi}_{4}$ must be replaced by the constraint that $u_{i,i}=0$. That the incompressibility condition should take the place of the pressure equation $\dot {p}=\delta H^{*}/\delta {\rm \pi}_{4}$ is consistent with the well-known result that the pressure usually serves as the Lagrange multiplier for the incompressibility constraint (Lanczos 1970; Badin & Crisciani 2018).

4.1. Hamilton–Jacobi equation

One of the most significant aspects of the Hamiltonian formalism is that it leads to the transformation theory of Hamilton (1833, 1834, 1835) and Jacobi (1837, 1842–1843) (see Whittaker (1904) and Lanczos (1970)), celebrated both for unifying particle mechanics with wave optics (Hamilton 1833) and for its relationship to the Schrödinger equation of quantum mechanics (Schrödinger 1926a,b). Here, we will obtain a Hamilton–Jacobi equation representing the Navier–Stokes problem.

In the context of discrete mechanics, Hamilton's principal function is obtained as the solution to the Hamilton–Jacobi equation, which is in turn defined by the functional form of the Hamiltonian. Hamilton's principal function provides the generating function for a canonical transformation to a new set of generalized coordinates and conjugate momenta for which the Hamiltonian vanishes identically, in which case Hamilton's equations do, in fact, become trivial. The new coordinates and their conjugate momenta are simply equal to their initial values.

In the present context, the role of Hamilton's principal function is played by a characteristic functional ${S}^{*}={S}^{*}[u_{i},p,t]$ (not to be confused with the action $\mathcal {S}^{*}$, although they are related; see Appendix A), which is the solution to the following Hamilton–Jacobi equation:

(4.10)\begin{equation} {H}^{*}\left[u_{i},p,\frac{\delta {S}^{*}}{\delta u_{j}},\frac{\delta {S}^{*}}{\delta p};t\right]+\frac{\partial {S}^{*}}{\partial t}=0, \end{equation}

where $\delta {S}^{*}/\delta u_{i}$ and $\delta {S}^{*}/\delta p$ are the Volterra (Volterra 1930) functional derivatives of ${S}^{*}$ with respect to the field quantities. Interested readers will find the derivation of (4.10) in Appendix A. Henceforth, we will refer to ${S}^{*}$ as ‘Hamilton's principal functional’. Substituting for the conjugate momenta in (4.5), we obtain for the Hamilton–Jacobi equation

(4.11)\begin{align} &\int \text{d}^{3}x\left[\frac{1}{2}\frac{1}{(\hat{\rho})^{2}}\frac{\delta {S}^{*}}{\delta u_{i}}\frac{\delta {S}^{*}}{\delta u_{i}}-\frac{1}{\hat{\rho}}\left(\hat{\rho}u_{i,j}u_{j} +{ p_{,i} - \mu u_{i,jj} - (\mu+\lambda)u_{j,ji}} - \hat{\rho}b_{i}\right)\frac{\delta {S}^{*}}{\delta u_{i}}\right. \nonumber\\ &\quad \left.+\,\frac{1}{2}\frac{1}{(\hat{\rho}')^{2}}\frac{\delta {S}^{*}}{\delta p}\frac{\delta {S}^{*}}{\delta p}-\frac{1}{\hat{\rho}'}\left({\hat{\rho}'p_{,i}u_{i} + \hat{\rho}u_{i,i}}\right)\frac{\delta {S}^{*}}{\delta p}\right]+\frac{\partial {S}^{*}}{\partial t}=0, \quad \hat{\rho}'\neq0. \end{align}

In contrast to the four original field equations – (3.1) and (3.2) – the Hamilton–Jacobi equation (4.11) is a single equation in Hamilton's principal functional ${S}^{*}$. This constitutes an equivalent formulation of the problem, as a complete and non-trivial solution to (4.11) is tantamount to an integration of Hamilton's equations (4.6a,b) and (4.7a,b) (note that it is not appropriate to take ${S}^{*}\equiv 0$ for the same reason that it is not appropriate to take $H^{*}\equiv 0$). In this way, we have reduced the problem of finding four separate field quantities to that of finding a single functional in those field quantities. One need only deduce (or even guess) the general form of ${S}^{*}$ in order to solve the problem. If an analytical expression for ${S}^{*}$ can be obtained, it will lead via canonical transformation to a new set of fields ($\phi _{i}, \phi _{4}$) and conjugate momenta which are simply equal to their initial values, giving analytical expressions for the four original fields $(u_{i},p)$.

The case of incompressible flow requires care, as $\hat {\rho }'\equiv 0$ and $\hat {\rho }'$ appears in the denominators of terms in (4.11). Even so, the Hamiltonian formulation remains well posed in the incompressible limit. Recall that, with $\hat {\rho }'\equiv 0$, $\mathcal {R}_{4}$ reduces to $\hat {\rho } u_{i,i}$, ${\rm \pi} _{4}$ vanishes identically and the Hamiltonian density reduces to (4.9). Substituting for the conjugate momenta ${\rm \pi} _{i}$ in (4.9), the corresponding Hamilton–Jacobi equation is

(4.12)\begin{align} \int \text{d}^{3}x\left[\frac{1}{2}\frac{1}{{\rho}^{2}}\frac{\delta {S}^{*}}{\delta u_{i}}\frac{\delta {S}^{*}}{\delta u_{i}}-\frac{1}{\rho}({\rho}u_{i,j}u_{j} + p_{,i} - \mu u_{i,jj} - {\rho}b_{i})\frac{\delta {S}^{*}}{\delta u_{i}}\right]+\frac{\partial {S}^{*}}{\partial t}=0, \quad \hat{\rho}'\equiv0, \end{align}

with $\delta {S}^{*}/\delta p=0$, since again ${\rm \pi} _{4}$ vanishes identically for incompressible flow. Here, the merit of starting from the compressible form of the equations becomes fully evident, as it would not necessarily have been clear that $\delta {S}^{*}/\delta p$ should vanish in the incompressible limit without knowing that in general ${\rm \pi} _{4}=\hat {\rho }'\mathcal {R}_{4}$. This is the form of the Hamilton–Jacobi equation as it relates to the traditional Navier–Stokes problem. In this case, the pressure is determined last of all, and is whatever it needs to be to enforce the incompressibility constraint $u_{i,i}=0$ (again consistent with the role of pressure as Lagrange multiplier Lanczos 1970; Badin & Crisciani 2018).

It must be acknowledged that the Hamilton–Jacobi equation developed above (either (4.11) for the compressible case or (4.12) in the incompressible limit) contains Volterra (1930) functional derivatives and is thus by no means trivial to solve. Indeed, it appears that solving such equations is itself a long-standing open problem in mathematics, which has received very little attention since the first half of the twentieth century (Michal 1926; Jordan & Pauli 1928; Levy 1951; Tatarskii 1961; Syavavko & Mel'nichak 1974; Dieudonne 1981; Koval'chik 1993). Nevertheless, if a rigorous theory of such equations can be developed, the present formulation of the Navier–Stokes problem might be solved as one special case. The present authors submit that such an endeavour is worthwhile and merits further study.

We conclude this section by remarking that, in the inviscid limit ($\mu =\lambda =0$), all of the preceding formalism remains perfectly well posed. In that limit, the present approach yields a mathematically equivalent second-order formulation of the inviscid Euler equations, as one would expect. Interested readers will find the full details in Appendix B.

5. Discussion

In this section we provide some qualitative interpretations of the developments of § 3. More specifically, we investigate the incompressible form (via constant, uniform density) of the Euler–Lagrange equations (3.9) and (3.10) when the residuals $\mathcal {R}_i$ and $\mathcal {R}_4$ are substituted.

Our motivation is again the simple example of § 1.1 for which the first-order non-Hamiltonian system $\dot {v}=-v$ was converted to the second-order Hamiltonian system $\ddot {v}=v$ by (manual) elimination of the non-conservative ‘damping’ term $\dot {v}$ (see Sanders (2022) for a similar result for the damped harmonic oscillator converting from second-order to fourth-order dynamics). Sanders (2023b) showed that the elimination process is ‘automated’ by the definition of the action in the first integral of (1.3), which is generalized to the action in (3.3) for our current continuum dynamics problem containing fields.

First consider the pressure equation (3.10) and corresponding natural boundary condition (3.15), which take the following incompressible forms:

(5.1)\begin{equation} {-}\mathcal{R}_{i,i}=0 \quad \forall x_j\in\mathcal{V}, \quad\text{subject to}\quad \mathcal{R}_in_i = 0 \quad \forall x_j\in\partial\mathcal{V}. \end{equation}

This higher-order field equation is simply the divergence of the residual $\mathcal {R}_i$. Upon substituting for $\mathcal {R}_i$ from (3.1) and subsequently imposing the incompressible continuity condition $u_{i,i}=\mathcal {R}_4/\rho =0$ from (3.2), we obtain

(5.2)\begin{equation} p_{,ii} ={-}[\rho u_ju_{i,j}]_{,i} + \rho b_{i,i}, \end{equation}

which is a Poisson equation for the pressure. The boundary condition is a Neumann type requiring the specification of the normal pressure gradient, $n_ip_{,i} = p_{,n} \equiv f(x_j,t)$, on the boundary, where

(5.3)\begin{equation} f(x_j,t) ={-}n_i[\rho\dot{u}_i+\rho u_ju_{i,j}-\mu u_{i,jj}-\rho b_i]. \end{equation}

Equation (5.2) and boundary condition (5.3) evolve the pressure in a manner that ensures the velocity field is solenoidal. This is a well-known pressure–velocity-based formulation commonly used in the numerical solution of incompressible flows (e.g. Ferziger & Peric 2002; Pozrikidis 2009).

Next, we consider the velocity equations (3.9) which, at present, have a more elusive physical interpretation. Here, we instead begin with the natural boundary conditions (3.13) and (3.14), which are due to the $\delta u_i$ and $\delta u_{i,j}$ variations. The incompressible versions of these equations are

(5.4a,b)\begin{equation} \rho \mathcal{R}_{i}u_{j}n_{j}+\mu \mathcal{R}_{i,j}n_{j}=0 \quad\text{and}\quad - \mu \mathcal{R}_{i}n_{j}=0 \quad \forall x_j\in\partial\mathcal{V}. \end{equation}

The boundary conditions involving the residual $\mathcal {R}_i$ are those compatible with the first-order Navier–Stokes equations, such as the no-slip and no-penetration conditions. Indeed, if we specify the velocity vector of the actual motion on the boundary, then $\mathcal {R}_i\equiv 0$ there. Note that the pressure of the actual motion on the boundary will be known from the simultaneous solution of (5.2).

However, the gradient terms $\mathcal {R}_{i,j}$ will introduce up to third-order spatial derivatives that must be specified. These represent the additional boundary conditions that must accompany the higher-order governing equation, which will be seen shortly to be second order in time and fourth order in space. Again, recall the example of § 1, whereby the system (1.2a–c) must be appended with a second (initial) condition specifying the (time) derivative of the coordinate $v(t)$. In the present context, these boundary conditions are ostensibly tantamount to specification of the viscous stress on the boundary by way of velocity gradients.

In general, the conditions at a boundary require two transition relations (Batchelor 2000; White 2006) to ultimately describe the momentum transport. Mathematically speaking, these conditions are the jump in velocity (momentum intensity) and the jump in stress (momentum flux). Under ordinary physical circumstances the velocity and stress are assumed to be continuous. However, this is one particular form of the transition relations, and there are familiar examples to which they do not apply. For example, at a liquid–gas interface the stress relation is modified to account for a non-zero jump in the normal stress that is balanced by a force due to surface tension (the tangential stress component usually still taken to be continuous). Similarly, in the event that molecular slip occurs, the typical transition relation gives an expression for the slip velocity (e.g. Thompson & Troian 1997; Thalakkottor & Mohseni 2016). In the case of energy transport, analogous conditions are needed regarding jumps in temperature (intensity) and heat flow (flux), which are recognized as the concept of thermal contact resistance.

We now turn our attention to the Euler–Lagrange equations (3.9), which upon imposing incompressibility and expanding derivatives of product terms yields

(5.5)\begin{equation} \rho\dot{\mathcal{R}}_i + \rho u_j\mathcal{R}_{i,j} = \rho \mathcal{R}_{j}u_{j,i}- \mu \mathcal{R}_{i,jj} \quad \forall x_j\in\mathcal{V}. \end{equation}

The left-hand side is the material derivative of the residual $\mathcal {R}_i$. Our purpose here is to observe which terms from the first-order Navier–Stokes equation are ‘eliminated’ in the higher-order formulation. Specifically, we are interested in the non-conservative viscous terms; while the body force $b_i$ could also be non-conservative, we will not concern ourselves with this possibility. Direct substitution of $\mathcal {R}_i$ into (5.5) generates many terms, but it is found that only one is cancelled: the viscous Laplacian of the (time derivative of the) velocity, namely $\mu \dot {u}_{i,jj}$. This term mutually appears from the $\rho \dot {\mathcal {R}}_i$ and $-\mu \mathcal {R}_{i,jj}$ terms in (5.5). To maintain notional clarity, we write the residual as

(5.6)\begin{equation} \mathcal{R}_i = \rho\dot{u}_i - \mu u_{i,kk} + \tilde{\mathcal{R}}_i, \end{equation}

where index $k$ has been used to avoid confusion with gradient operators in (5.5) having index $j$, and $\tilde {\mathcal {R}}_i = \rho u_ku_{i,k}-\rho b_i$ are the remaining terms in the residual. Substituting the above into the first and last terms of (5.5), cancelling the aforementioned $\mu \dot {u}_{i,jj}$ term and then dividing out by the density gives

(5.7)\begin{equation} \rho\ddot{u}_i+\dot{\tilde{\mathcal{R}}}_i+u_j\mathcal{R}_{i,j} = \mathcal{R}_ju_{j,i} - \nu[-\mu u_{i,kkjj} + \tilde{\mathcal{R}}_{i,jj} ], \end{equation}

where $\nu =\mu /\rho$ is the kinematic viscosity (recall that all variables are non-dimensional). We see that this equation is second order in time and fourth order in space. Viscous terms still appear in the equation including second- and third-order spatial derivatives. Nevertheless, the technique detailed by Sanders (2023b) and employed here evidently ensures that (5.7) has a corresponding Hamiltonian structure.

6. Case study

We can explore how this method can be applied by considering a simplified example with a known field solution. In looking at the variety of cases in which the Navier–Stokes equations have a known analytical solution, the simplest are those involving steady flows. While the Euler–Lagrange equations (3.9), (3.10) can be written for these cases, the corresponding Hamilton–Jacobi equation is trivial because for steady flows the fields are already equal to their initial values.

It is therefore worthwhile to examine the simplest unsteady flows, which should result in a non-trivial Hamilton–Jacobi equation. Indeed, there exists a class of flows for which the Navier–Stokes equations take the same simplified form: those in which the flow is incompressible and unidirectional (Batchelor 2000). This class of problems include both of Stokes's flows (Stokes 1851), in which a semi-infinite fluid is influenced by a boundary moving in its own plane. In the first of these cases, the boundary is impulsively started and in the second, the boundary oscillates. We can also include developing flow in a channel or pipe. The only difference between these flows results from initial and boundary conditions, but the Navier–Stokes equations and therefore the present Hamilton–Jacobi equation take the same form.

Here we will examine the case in which there is motion only in the $x_1$ direction, and the velocities take the form $\{u_{i}\} = \{u_{1}(x_{2},t),0,0\}$. In the absence of a body force, our pressure gradient in the $x_1$ direction is solely a function of time and the pressure gradients in the $x_2$ and $x_3$ directions are zero. There are thus only two unknown field quantities: $u_{1}(x_{2},t)$ and $p(x_{1},t)$, where $p$ is linear in $x_{1}$. The field equation of primary interest is

(6.1)\begin{equation} \mathcal{R}_{1} \equiv \rho \dot{u}_{1} + p_{,1} - \mu u_{1,22} = 0, \end{equation}

and the remaining field equations are satisfied automatically by the assumed form of the fields. Following the procedure described above, the momenta conjugate to $u_{1}$ and $p$ are given by

(6.2a,b)\begin{equation} {\rm \pi}_{1}\equiv \rho\mathcal{R}_{1}, \quad {\rm \pi}_{4}\equiv0.\end{equation}

This results in a Hamiltonian density given by

(6.3)\begin{equation} \mathcal{H}^{*}=\frac{1}{2}\frac{1}{{\rho}^{2}}{\rm \pi}_{1}{\rm \pi}_{1}-\frac{1}{\rho}(p_{,1}-\mu u_{1,22}){\rm \pi}_{1}. \end{equation}

Hamilton's principal functional ${S}^{*}={S}^{*}[u_{1},p,t]$ can be expressed as an integral over $x_{2}$ only, since the other spatial coordinates do not appear and may be integrated out. In this way, we may write the Hamilton–Jacobi equation as follows:

(6.4)\begin{equation} \int \text{d}\kern 0.06em x_2\left[\frac{1}{2}\frac{1}{{\rho}^{2}}\frac{\delta {S}^{*}}{\delta u_{1}}\frac{\delta {S}^{*}}{\delta u_{1}}-\frac{1}{\rho}(p_{,1}-\mu u_{1,22})\frac{\delta {S}^{*}}{\delta u_{1}}\right]+\frac{\partial {S}^{*}}{\partial t}=0, \end{equation}

with $\delta {S}^{*}/\delta p=0$. The solution to (6.4) would provide a canonical transformation to a new set of coordinates, giving analytical expressions for $(u_{1},p)$.

Despite knowing the analytical solution for these fields in this particular example, the present authors have not been able to solve this Hamilton–Jacobi equation, since again the solution of such equations is itself an open problem (Michal 1926; Jordan & Pauli 1928; Levy 1951; Tatarskii 1961; Syavavko & Mel'nichak 1974; Dieudonne 1981; Koval'chik 1993). This example therefore appears to be a good place to start for tackling the general problem. Another interesting example to consider might be a two-dimensional Taylor–Green vortex such as that considered by Wu, Ma & Zhou (2006).

7. Conclusion

This paper has presented a novel Hamiltonian formulation of the isotropic Navier–Stokes problem for both compressible and incompressible fluids. This canonical formulation opens several previously unexplored avenues toward a final resolution of the problem, which we briefly describe below.

Perhaps the most obvious route would be to solve the Hamilton–Jacobi equation – either (4.11) for the compressible case or (4.12) for the incompressible case – for Hamilton's principal functional ${S}^{*}[u_{i},p,t]$ directly. If a complete solution for ${S}^{*}$ can be found, it will lead via canonical transformation to a new set of fields which are equal to their initial values, thereby giving analytical expressions for the original velocity and pressure fields. Alternatively, if one can simply establish based on emerging analytical techniques that a complete solution to this Hamilton–Jacobi equation does (or does not) exist under the usual assumptions, that will also settle the question of existence of solutions.

An alternative strategy might be to investigate the corresponding Lagrangian formulation based on the action $\mathcal {S}^{*}$ as given by (3.3). Because the first- and second-order formulations are mathematically equivalent (recall the proof in § 3.2), $\mathcal {S}^{*}$ must have as many local minima as there are solutions to the traditional, first-order formulation. Intuitively, it seems as though it ought to be possible to determine – or at least to establish bounds on – the number of critical points an action has based on the form of the Lagrangian (Van den Berg & VanderVorst 2002; Kalies & VanderVorst 2004). If one can establish that, under the usual assumptions, $\mathcal {S}^{*}$ always has exactly one local minimum, or else demonstrate that there are cases where it fails to achieve a local minimum, that too will resolve the question of existence and uniqueness.

By no means is either of the above programs trivial. As pointed out in § 4.1, solving equations containing Volterra (1930) functional derivatives is itself a long-standing open problem in mathematics, which has received very little attention since the first half of the twentieth century (Michal 1926; Jordan & Pauli 1928; Levy 1951; Tatarskii 1961; Syavavko & Mel'nichak 1974; Dieudonne 1981; Koval'chik 1993). One might even go so far as to call it a ‘forgotten’ open problem (as did one of the reviewers of the present paper, who generously drew our attention to Jordan & Pauli 1928; Levy 1951; Tatarskii 1961; Syavavko & Mel'nichak 1974; Dieudonne 1981; Koval'chik 1993). We see the lack of work on such equations as a challenge, yes, but at the same time we also see it as a significant opportunity for advancing the field of analytical continuum mechanics. Perhaps, despite an apparent increase in complexity, a rigorous theory of such equations can be developed after all, in which case the present formulation of the Navier–Stokes problem might be solved as one example. We submit that, at the very least, such an endeavour merits further study, which we intend to continue in future work.

Finally, it is worth noting that the techniques employed here are by no means specific to the Navier–Stokes problem, nor are they restricted to the field of classical mechanics. The suitably averaged principle of least squares (Sanders 2021, 2022, 2023a,b; Sanders & Inman 2023) can be applied to any traditionally non-Hamiltonian dynamical system in order to formulate a mathematically equivalent higher-order Hamiltonian system. It is believed that this fundamental result will also find uses in other branches of pure and applied mathematics.