Drawing an analogy with quantum mechanics, a new Lagrangian is proposed for a variational formulation of the Navier–Stokes equations which to-date has remained elusive. A key feature is that the resulting Lagrangian is discontinuous in nature, posing additional challenges apropos the mathematical treatment of the related variational problem, all of which are resolvable. In addition to extending Lagrange's formalism to problems involving discontinuous behaviour, it is demonstrated that the associated equations of motion can self-consistently be interpreted within the framework of thermodynamics beyond local equilibrium, with the limiting case recovering the classical Navier–Stokes equations. Perspectives for applying the new formalism to discontinuous physical phenomena such as phase and grain boundaries, shock waves and flame fronts are provided.
For physical systems formulated within the framework of Lagrange's formalism, the dynamics are completely defined by only one function: the Lagrangian. This methodical concept successfully applies to, for example, conservative Newtonian mechanics. Contrary to this, in continuum theories many open problems remain unsolved to this day, especially when considering dissipative processes; the viscous flow of a fluid, given by the Navier–Stokes equations, is a typical example where a formulation by a Lagrangian is missing.
Many attempts at finding a variational formulation of the Navier–Stokes equations have been made: Millikan  performed an investigation by assuming a Lagrangian of the form in terms of the velocity u, the pressure p and their first-order derivatives. Despite his rigorous treatment of this inverse problem, his result applies to special flow geometries only. Consequently, in general a different approach is required based on the representation of the observable fields by potentials, i.e. by auxiliary fields representing the observables. For the particular case of inviscid flows, Clebsch  was successful in this by means of the potential representation 1.1 known as the Clebsch transformation [3,4,5]. There are various field theories, for instance Maxwell theory, where a representation of observables in terms of potentials is required in order to establish a proper variational formulation. At first glance, it seems to be a matter of experience that potential fields are required for finding Lagrangians in different field theories. However, in Scholle  a convincing explanation is given as to why in continuum theories the use of potentials is absolutely necessary for the construction of a Lagrangian: in order to fulfill the invariance with respect to the full Galilean group, at least one field must be non-measurable and therefore be a potential. In the same paper, a general scheme for Lagrangians is constructed. Using Noether's theorem, canonical formulae give rise to the identification of the relevant observable fields like mass density and flux density, momentum density, stress tensor, energy density and Poynting vector.
It should also be noted that the existence of the Clebsch variables ζ,α,β for an arbitrary velocity field u can be taken for granted only locally. Their global existence depends on the topological features of the flow; for details we refer, for example, to [7,8]. On the other hand, completeness of the Clebsch representation may be reached by additional pairs of variables. Here we restrict our theory to the classical form (1.1).
Since viscosity leads to dissipation and therefore to the irreversible transfer of mechanical energy to heat, thermal degrees of freedom have to be considered in order to remain consistent to Noether's theorem which implies conservation of energy for systems with time-translation invariance, because otherwise the time-translation invariance would have to be violated by an explicit time-dependence. Seliger & Whitham  made a suggestion as to how to embed thermal degrees of freedom in a variational formulation of fluid flow via the Lagrangian 1.2 and 1.3 which depends on the specific inner energy e(ϱ,s), given in terms of the mass density ϱ and the specific entropy s, the three Clebsch potentials ζ, α, β and an additional potential field ϑ. The meaning of the latter becomes apparent by calculating the Euler–Lagrange equation with respect to s, giving the ‘potential representation’ 1.4 for the temperature T, used three decades previously by Van Dantzig, who termed the field ϑ as thermasy . Although still restricted to adiabatic and therefore reversible processes, the Lagrangian (1.2) represents a momentous step forward because of the rudimentary embedding of thermodynamics.
By comparing the potential representation (1.3) with the one proposed by Clebsch (1.1) for the isothermal case, it becomes apparent that any kind of extension of the system, by additional degrees of freedom as well as by additional physical effects, requires an adjustment of the potential representation; see, for example, Wagner . Scholle  provides a general explanation for the necessary use of different potential representations of the observables for different physical systems along the lines of a rigorous analysis of the fundamental symmetries the Lagrangian has to fulfil, with particular regard to Galilean invariance. In the same paper, an easily manageable symmetry criterion for verifying Galilean invariance is derived. Based on these preliminary findings, Scholle [12,13] suggested a Lagrangian for viscous flow by supplementing the Lagrangian (1.2) with additional terms leading to partial success: on the one hand, the phenomenon ‘viscosity’ occurs in a qualitatively correct manner, as demonstrated by three simple flow examples. On the other hand, the equations of motion resulting from the variation of Hamilton's principle differ from the Navier–Stokes equations and, therefore, their solutions reveal notable quantitative differences to those of the latter. Similar experiences are reported by Zuckerwar & Ash , who made an analogous suggestion for a Lagrangian considering volume viscosity in particular.
Despite this partial success, the need to improve the existing approach is obvious in order to obtain solutions from Hamilton's principle suitable for relevant flow problems. In this article, an innovative idea by Anthony  is applied, based on the reformulation of the Lagrangian in terms of complex fields. This can also be understood as the inversion of Madelung's idea  of reformulating the complex Schrödinger's equation into a hydrodynamic form.
2. Construction of the Lagrangian
First, as demonstrated in prior work [17,18,12,13], Seliger and Whitham's Lagrangian (1.2) can be re-written alternatively as 2.1 in terms of the extended set of independent fields ψ=(u,ζ,α,β,ϱ,s,ϑ) and their material time derivatives 2.2 The above form of the Lagrangian yields two benefits: first, the potential representation (1.3) of the velocity field results from a variation with respect to u and hence does not need to be prescribed; second, by adding terms to the Lagrangian depending on first-order derivatives of u in order to consider friction, the extended Lagrangian remains of first order. The latter is a useful feature because otherwise, i.e. in the case of a Lagrangian containing second-order derivatives of the fields, the computation of (i) the corresponding Euler–Lagrange equations and (ii) the canonical densities and flux densities resulting from Noether's theorem become more complicated. It is also very useful to avoid derivatives of order higher than one when applying Ritz's direct method to problems formulated in curvilinear coordinates.
2.1 Conventional approach and examples
For the extension of the Lagrangian (2.1) to incorporate viscosity, it seems reasonable to simply add terms modifying the entropy balance resulting from variation with respect to the thermasy ϑ as 2.3 the homogeneity of this expression indicating that only adiabatic processes are considered in (2.1). Note that the above balance alternatively results from Noether's theorem with respect to the transformation ϑ→ϑ+ε with ε=const., which is a symmetry transformation of the Lagrangian (2.1). In order to take into account the production of entropy, the symmetry with respect to the transformation ϑ→ϑ+ε has to be broken, which is achieved in the easiest way by adding a term linearly dependent on ϑ, i.e. ϑϕd/T, to the Lagrangian, where the dissipation heat ϕd should be positive and depend on the spatial derivatives of the velocity as the primary cause for the physical phenomenon ‘viscosity’. Both are satisfied by assuming for ϕd a quadratic dependence on ∂jui, according to the classic literature on viscous flow [20,21]. Finally, via the factor 1/T the character of the entropy as ‘weighted heat’, according to δQ=TdS, is accounted for. The above considerations provide the motivation behind the following extended Lagrangian [12,13]: 2.4 where η is the shear viscosity, η′ the volume viscosity of the fluid and 2.5 is the tensor of the shear rate; tr denotes the trace of a tensor. The temperature T, according to classical thermodynamics, is given by 2.6
Note that neither an external force nor heat conduction is considered here. Now, by variation with respect to the thermasy ϑ, the equation 2.7 results as an entropy balance with an entropy production rate on the right-hand side due to dissipation, as expected. Furthermore, the above Lagrangian fulfils the necessary symmetry requirements for Galilean invariance, as analysed in detail in §A.1 of appendix A; however, an unexpected feature transpires. The momentum density p, resulting as a canonical Noether observable, does not equal the mass flux density ϱu. The difference between both 2.8 termed quasi-momentum density, needs to be explained physically. According to Scholle , p* could be due to contributions to the system's momentum balance beyond the scope of the continuum hypothesis on a molecular scale, e.g. Brownian motion. This question is discussed in more detail in §5. Regardless, the dynamics induced by the Lagrangian (2.4) go beyond the scope of classical theory, the resulting equations of motion differing significantly from the Navier–Stokes equations. In the case of incompressible flow and negligible buoyancy they read [12,13,19] 2.9 2.10 2.11 with constant mass density ϱ=ϱ0 and kinematic viscosity ν:=η/ϱ0. Looking at the above set of PDEs, two striking features immediately become apparent, namely:
(i) the resulting field equations are third-order PDEs, not second-order ones like the Navier–Stokes equations;
(ii) the thermasy ϑ appears explicitly in the field equations as an additional degree of freedom.
The above qualitative features have also been found in the case of compressible flow with pure volume viscosity by Zuckerwar & Ash , and the appearance of an additional physically relevant degree of freedom, in particular, appears also in the variational formulation of heat conduction proposed by Anthony  and is interpreted by him as a measure for the deviation from the thermodynamical local equilibrium. Similar assumptions are made in Zuckerwar & Ash . At first glance, there seems to be a realistic chance to interpret the additional terms and degrees of freedom in the above evolution equations (2.7), (2.9)–(2.11) as an extension of the classical theory towards non-equilibrium thermodynamics. In order to test this hypothesis, three ‘benchmark tests’ have been performed [12,13]:
(i) Plane Couette flow, i.e. shear flow between two parallel plates, the upper one of which is moving with a constant speed U, representing one of the most prominent examples in fluid dynamics. Considering the unidirectional flow u=u(y)ex with boundary conditions u(0)=0 and u(h)=U, a solution of equations (2.9)–(2.11) is given by the linear velocity profile u(y)=Uy/h and a constant pressure p in full accordance with the Navier–Stokes equations.
(ii) As an example of a transient flow, a suddenly moving plate is considered: a horizontal plate of infinite extension is covered by a fluid at rest at t<0, see figure 1a. At t=0 the plate suddenly starts moving with constant velocity U in the horizontal direction, invoking a flow inside the fluid. The initial conditions, t0=0 and u0=0, have to be considered. Since no characteristic length is contained in this problem, a representation of the velocity profile in terms of a similarity variable u=Uf(ξ) with  is used, leading to the solution of the equations of motion (2.9)–(2.11), whereas the classical solution of the Navier–Stokes equations reads f(ξ)=1−erf(ξ/2)  involving the error function erf. In figure 1a, both resulting profiles are compared to each other, revealing qualitatively concordant profiles with quantitative differences.
(iii) The Lamb–Oseen vortex , the flow geometry of which is shown in figure 1b, is another example of a transient flow: using cylindrical coordinates r,φ,z, a similarity variable is available allowing for a representation of the velocity as u=u(r,t)eφ with the time-depending profile u(r,t)=Γf(ξ)/(2πr) for a given circulation Γ. According to Scholle [12,13,19], from (2.9)–(2.11) the analytical form f(ξ)=1−ξK1(ξ) is obtained with the modified Bessel function of first order, K1. In contrast, occurs in the classical solution . Both solutions are compared in figure 1b, revealing again qualitative accordance with quantitative differences.
Note that for the above three examples (i)–(iii), ϑ/T=t has been used as a particular solution of the evolution equation (2.11), leading to the simplified form 2.12 of the equations of motion (2.9). Summarizing the above examples, it can be ascertained that the phenomenon of viscosity is at least qualitatively captured by Hamilton's principle. Hence, the Lagrangian (2.4) reflects a relevant step forward to a satisfying description of viscous flow within the framework of Lagrange formalism. On the one hand, the differences between the classical theory and the results obtained by variation are quite pronounced for the transient flow examples (ii) and (iii). It is open to dispute if they can be explained as non-equilibrium effects. On the other hand, alternative equations of motion containing additional terms related to thermal relaxation are discussed in the literature; for instance, we refer to Ash & Zardadkhan [22,23] who made use of the Navier–Stokes equations supplemented with a ‘pressure relaxation term’, implying a vortex solution considerably different from the classical Lamb–Oseen vortex in order to resolve discrepancies between existing models and observations on dust devils and tornados. Apart from the examples reported in Scholle [12,13,19], we consider three further examples in order to get things straight concerning the question of whether the non-classical form of the equations of motion (2.9)–(2.11) can consistently be explained as being due to non-equilibrium thermodynamics. These are:
(iv) The Taylor–Couette flow between two cylinders of radius ri and ra, see figure 2a, invoked by a rotation of the inner cylinder with angular velocity ω0, is characterized by a flow geometry u=u(r)eφ. By applying this to the equations of motion (2.10) and (2.12) and considering the boundary conditions u(ri)=ω0ri and u(ra)=0, we obtain as a solution and in perfect agreement with classical theory.
(v) Plane Poiseuille flow between two parallel plates, a distance h apart, driven according to figure 2b by a pressure gradient (p1−p2)/l, with a unidirectional flow geometry u=u(y)ex assumed. Considering no-slip conditions u(0)=u(h)=0 at the lower and upper plate, the solution of the equations of motion (2.10), (2.12) reads 2.13 and 2.14 By identifying K=(p1−p2)/(ηl), the above velocity profile (2.13) is in perfect agreement with the classical solution . However, the associated pressure (2.14) contains an additional term ηKu(y)t, by which the adherence of the boundary conditions for the pressure p at the inflow and the outflow is inhibited. Moreover, the pressure is unsteady and, as a non-physical feature, it tends to infinity with increasing time. The reason for the latter problem stems from the choice of the particular thermasy solution ϑ/T=t of the evolution equation (2.11) which is increasing with time as well, a fact that seems to be problematic not only for this specific flow problem but for problems in fluid mechanics in general, as stated in Scholle . In their response to the comment , Zuckerwar and Ash  suggested to construct a time-independent solution of the evolution equation (2.11), fulfilling u⋅∇(ϑ/T)=1. Following their suggestion, as a steady solution for the weighted thermasy the expression 2.15 is obtained with arbitrary integration function f1(y), see appendix A.2. The associated solution of the equations of motion (2.9) and (2.10) reads, according to appendix A.2, implicitly as 2.16 2.17 2.18 with the Gaussian hypergeometric function 2F1. The velocity profile (2.17) is visualized in figure 2b (right) and differs markedly from the classical parabolic profile (2.13) (left), especially due to the fact that its first-order derivative vanishes at the walls, u′(0)=u′(h)=0, thus indicating a zero wall shear stress. As another conspicuous feature, the y-dependence of the pressure inhibits again fulfilment of the boundary conditions for the pressure p at the inflow and the outflow. In summary, no physically meaningful steady solution of (2.9)–(2.11) can be constructed for a plane Poiseuille flow.
(vi) Poiseuille flow in a curved channel, see figure 2c, as investigated by Richter , who discovered qualitatively the same problems known from the plane Poiseuille flow. In particular, the resulting pressure solution, is found to similarly contain a non-physical term increasing with time. Other attempts to find a time-independent solution matching the boundary conditions have remained unsuccessful, as in the case of the previous example.
Summarizing the above benchmark tests, we conclude that only for two of the six examples, namely the shear-driven flows (i) and (iv), are the classical solutions recovered. For the two transient flows, (ii) and (iii), reasonable solutions of the field equations and boundary conditions have been found, with quantitatively different velocity profiles compared to the classical solutions. For the two pressure-driven flows, (v) and (vi), no adequate solutions of the field equations could be constructed which simultaneously fulfil the pressure boundary conditions. Hence, the variational principle based on the Lagrangian (2.4) does not recover the dynamics of viscous flow in a proper way, since its applicability seems to be restricted to special flow problems only.
Moreover, four of the six benchmark solutions contradict the hypothesis that the differences between them and classical theory can be explained by effects beyond the scope of thermodynamic equilibrium: apart from the non-physical features discovered above, one would expect that the non-equilibrium solution tends towards the classical equilibrium solution if a special relaxation parameter in the problem, physically related to the deviation from equilibrium, tends to infinity. This is not the case here since no additional parameters exist but mass density, viscosity and specific heat.
Analysing the above examples in more detail, the explicit appearance of the weighted thermasy ϑ/T in the equations of motion (2.12) seems to be the crucial issue leading to non-physical solutions, since ϑ/T turns out to be of unlimited growth, either spatially or temporally, which also prohibits its interpretation in connection with non-equilibrium thermodynamics. Moreover, the anomalous relation (2.8) between mass flux density and momentum density implies that the discrepancy between mass flux and momentum also tends to increase spatially or temporally due to the ϑ/T-dependence.
A modification of the Lagrangian (2.4) is provided below which overcomes the above-highlighted anomalies.
2.2 Alternative approach based on complex fields
In 1927, Madelung  discovered a remarkable analogy between quantum mechanics and fluid mechanics by reformulating the complex Schrödinger's equation into a hydrodynamic form: by decomposing the quantum mechanical state function ψ into modulus and phase according to 2.19 Schrödinger's equation is transformed to the set of PDEs 2.20 and 2.21 These are obviously the equations of motion of a kind of fluid, the so-called Madelung fluid: equation (2.20) is the continuity equation and (2.21) is Bernoulli's equation for a fluid with the ‘unusual’ pressure function and with vorticity-free velocity field u=∇Φ. Based on these substitutions, Madelung established a fluid mechanics picture of Schrödinger's theory.
Many years later, Anthony  suggested the reverse of this idea, i.e. form a ‘Schrödinger-picture’ of fluid mechanics and thermodynamics, by combining the density ϱ and the Clebsch variable ζ in (1.2) to form a complex matter field ψ according to (2.19). Moreover, he introduced two more complex fields, namely a complex vortex potential Ω by combining the two remaining Clebsch variables α,β and the complex field of thermal excitation χ, giving the temperature by its absolute square: . The motivation for this transformation is originally given by Anthony's entropy concept: the entropy balance results from a canonical procedure related to the phase translation invariance of the complex fields as a balance of second kind within the framework of second variation and related stability criteria. Furthermore, Anthony states that by the complex representation a basic concept is given for an accurate formulation of thermodynamics of irreversible processes within the framework of Lagrange formalism.
For convenience, we only apply a partial transformation to complex fields in a slightly modified form to the Lagrangian (2.4): by introducing T0 as a constant reference temperature and c0 as a reference constant for the specific heat and considering the identity motivation is given for the generalized definition of the field of thermal excitation as 2.22 supplemented by the substitution 2.23
for the Clebsch variable ζ. Note that in (2.22) another constant, ω0, has been introduced due to dimensional reasons, like T0 and c0 before. Although there is no general rule how the three constants T0,c0 and ω0 have to be chosen, it is reasonable to choose the reference temperature T0 as a ‘typical’ temperature and c0 as a ‘typical’ specific heat. In the case of an incompressible flow with constant specific heat, discussed subsequently in §(c), the choice of c0 is obvious and the choice of T0 is arbitrary since the resulting Lagrangian (2.31) does not depend on it any more. By contrast, the choice of ω0 is not obvious, nor how the physics is affected by it. This will be analysed and discussed carefully in the following.
On substituting for (2.22) and (2.23), it follows that 2.24 and 2.25 with the sawtooth function (see also figure 3) 2.26 While the first equation (2.24) allows for a unique transformation of the respective real-valued terms in the Lagrangian (2.4) into terms of the complex thermal excitation field, the second equation (2.25) reveals that no equivalent for the thermasy ϑ explicitly appearing in the friction term of (2.4) can be constructed in terms of the complex field χ. The reason for this is the non-uniqueness of the argument of a complex number. The obviously most feasible way to resolve this issue is the use of as a substitute for ϑ, leading to the modified Lagrangian 2.27 Comparing this Lagrangian with (2.4), two remarkable differences are discernible: first, the modified Lagrangian is discontinuous due to the logarithmic term; second, the angular frequency ω0, which has primarily been introduced for dimensional reasons, becomes a relevant parameter, the physical meaning of which will be clarified subsequently. The most striking feature, however, is that the unlimited weighted thermasy ϑ/T appearing in (2.4) has been replaced by an expression with limited values between −π and π.
2.3 Incompressible flow with constant specific heat and external force
We consider a fluid flow with constant mass density, constant specific heat and without thermal expansion: 2.28 2.29 2.30 Note that (2.29) implies for the field of thermal excitation (2.22) in accordance with Anthony . Furthermore, since for the incompressible case volume viscosity is excluded from the very beginning, the term with η′ in (2.27) is absent. Finally, we add the specific potential energy V =V (x,t) of an external force, leading to the simplified Lagrangian: 2.31
3. Variation with discontinuous Lagrangian: general formalism
In conventional variational calculus, Euler–Lagrange equations can be computed if the Lagrangian is two times continuously differentiable . If this basic requirement is not fulfilled, a non-standard approach is required for variation, which is developed in the following.
We consider a variational principle δI=0 where I is given by 3.1 depending on independent fields ψi with i=1,…,N and ψN=φ. The Lagrangian ℓ is assumed to be discontinuous with respect to φ at fixed values φn, n=1,…,NS, but continuously differentiable with respect to all other fields, ψi with i<N, and also continuously differentiable with respect to all derivatives of any field, including and ∇φ. In three-dimensional space, the discontinuities with respect to φ become manifest along surfaces Sn(t) defined by 3.2 intersecting the system's volume V into a finite number NS+1 of sub-volumes according to: where the sub-volume V n denotes the region between Sn and Sn+1 apart from V 0 and V NS denoting the region between the system's boundary ∂V and S1 or SNS, respectively.
From the physical viewpoint, these time-dependent interfaces, Sn, may be related to any kind of discontinuous phenomena like phase boundaries between non-mixable fluids, propagating shock fronts in gaseous media or flame fronts. Their local propagation velocity is denoted by vs and the vector normal to the surface by n (see figure 4); the orientation of n is defined by the convention n⋅vs>0.
3.1 Euler–Lagrange equations
First, only the subset of variations δψi=0 is considered with δψi=0 at the interfaces Sn and at the system's boundary ∂V . Free variation is assumed inside the sub-volumes V n. Note that this kind of variation does not cause any shift of the interfaces Sn. Under these assumptions, the usual derivation procedure leading to the Euler–Lagrange equations can be performed separately inside each sub-volume, and hence the well-known Euler–Lagrange equations , 3.3 remain valid piecewise at each sub-volume V n.
3.2 Matching conditions
Next, a larger set of variations with δψi≠0 at the interfaces Sn is considered for i=1,…,N−1, but the variation of ψN=φ is again restricted to δφ=0 in order to exclude any shift of the interfaces themselves. Like before, δψi=0 is again prescribed at the system's boundary. Now, variation with respect to ψ1,…,ψN−1 leads to 3.4 using the abbreviation in (3.3) for the Euler–Lagrange expressions. By means of the Gaussian theorem and Reynold's transport theorem well-known from fluid dynamics  with vs being the velocity of the propagating interface Sn (figure 4), the variation takes the form 3.5 Since δψi=0 at ∂V n and at the initial and final time t1,2, the above identity simplifies, making use of the Euler–Lagrange equations (3.3), to In the following, the limit of the respective discontinuous expression by approaching it from the front side (subscript +) or the back side (subscript −) of the interface Sn is indicated by [⋯ ]±. Front and back side are defined by the orientation of the normal vector n according to figure 4. Then, by decomposing each surface integral over ∂V n in one integral along the front side of Sn and another one along the back side of Sn+1, we find in general 3.6 and, in particular, where the double square bracket indicates the jump at the interface: . Thus, variation δI=0 delivers 3.7 as matching conditions for the generalized fluxes, i=1,…,N−1 at each interface.
Independent of the formal proof given above the matching conditions can also be understood as natural boundary conditions at the phase boundaries in a multiphase flow when assuming that all phases of the flow consist of the same liquid, leading to the same equation (3.7).
3.3 Production condition
Finally, δψi=0 is again prescribed at the system's boundary, but apart from this free variation of all fields is allowed overall inside V , including free variation of ψN=φ. As a consequence of the latter, the position of the interfaces Sn is varied too, see figure 5: an arbitrary point x of Sn defined by (3.2) is shifted to a different position x+δx according to
leading to the identity: 3.8 for the local shift δx of the discontinuity interface Sn. In the case that Sn is shifted in the forward direction, n⋅δx>0, in a thin layer of thickness n⋅δx the value of ℓ is changed from [ℓ]+ to [ℓ]−, hence the integral (3.4) has to be supplemented by the following correction term: Since ∇φ∥n, the identity nδφ=−(δx⋅∇φ)n=−(δx⋅n)∇φ results, leading to Variation with respect to φ leads to the jump condition, 3.9 for the flux related to φ at the interface Sn. From the physical viewpoint, this is related to the production of the associated integral quantity. Therefore, we call (3.9) the production condition.
There is an alternative way to derive the above production condition (3.9) along the lines of distribution theory sketched in §A.5 of appendix A, showing that the generalized form of the formalism can be understood in terms of standard Lagrange formalism.
4. Resulting equations of motion and matching conditions
4.1 Equations of motion
In §A.3 of appendix A the Euler–Lagrange equations resulting from variation with respect to the elementary fields are calculated. Based on this, in §A.4 the corresponding equations of motion are derived. As a result, the equations of motion for the observable fields (A.15), (A.20) are 4.1 and 4.2 where fn.e. is used as an abbreviation for 4.3 Obviously, by (4.1) the continuity equation is perfectly reproduced, whereas the equations of motion (4.2) differ from the original Navier–Stokes equations by some additional forces fn.e.. Since according to (4.3) the latter contain the complex field of thermal excitation, χ, the corresponding evolution equation (A.19), 4.4 has to be considered additionally in order to complete the set of equations.
Reconsidering the aforementioned hypothesis of Anthony  as well as Zuckerwar & Ash  to identify the additional forces occurring in the equations of motion as contributions due to a deviation from the thermodynamic local equilibrium, one would expect a limit case leading to the classical dynamics, as already discussed at the end of §(a). Indeed, according to (4.3) the extra forces are scaled down when increasing the parameter ω0 and the physical dimension of ω0 is a reciprocal of time, suggesting its interpretation as a relaxation rate. This interpretation is underpinned by considering the limit , leading to vanishing of the extra forces, fn.e.→0, and therefore to a full reproduction of the Navier–Stokes equations by means of (4.2). In this equilibrium limit, the evolution equation (4.4) becomes meaningless, since the field of thermal excitation does not appear in the equations of motion any more.
It is noteworthy that the limit can be applied successfully to the equations of motion (4.2) on the one hand, reproducing Navier–Stokes equations, but cannot be applied directly to the Lagrangian (2.31) on the other hand, most likely because the mechanical equations (4.1) and (4.2) are decoupled from thermodynamics, whereas in the viscosity term of the Lagrangian (2.31) the occurring mechanical and thermodynamic degrees of freedom are strictly coupled.
For finite but sufficiently large values of ω0, the additional forces fn.e.→0 due to thermodynamic non-equilibrium remain small compared with viscous, external and pressure forces. According to (4.3), they consist of a factor expeditiously fluctuating between −π and π, of quadratic terms with respect to velocity gradients and of third-order derivatives of the velocity.
4.2 Matching conditions
As shown in §(b), variation with respect to the elementary fields, except for χ, induces matching conditions (3.7) at each interface. These are 4.5 4.6 4.7 4.8 According to the first condition (4.5), the normal component of mass flux density has to be continuous, which physically corresponds to the conservation of the mass passing the interface. By inserting (4.5) into condition (4.7), it reduces to [[α]]=0, implying continuity of the Clebsch variable α. In order to understand the physics behind condition (4.8), we have to take into account that at each discontinuity the phase of the thermal excitation jumps from π to −π or vice versa. In any case, the sign of changes when turning from the back side to the front side of an interface. As a consequence, condition (4.8) entails 4.9 which implies a reversal of the direction of the shear rate vector at the inner boundary. Physically, this is associated with a slip occurring at the interface and can again be interpreted as a phenomenon due to a deviation from thermodynamic equilibrium.
4.3 Production condition and thermodynamic aspects
In order to apply formula (3.9) for the production condition, the complex field of thermal excitation has to be decomposed into modulus and phase according to , leading to the real-valued form 4.10 of the Lagrangian (2.31), where S again denotes the sawtooth function. Then, the production condition reads: 4.11 Despite the reversal of the shear rate tensor at the interfaces according to (4.9), its square, , remains continuous. Hence it follows that , leading to: 4.12 The above condition reveals a discontinuity in the flux of the inner energy and therefore the production of inner energy due to dissipation at the interfaces. The latter one can alternatively be related to the volume by estimating the gradient of the thermal phase as ∇φ≈(2π/d)n, where d denotes the distance between two interfaces. As a consequence, the inhomogeneity at the right-hand side of the balance (4.12) can, according to , be re-interpreted as the mean production of inner energy related to the volume, at least in the sense of a statistical treatment.
Another source of inner energy production is already given by equation (A.21): 4.13 which takes the form of a classical balance equation in continuum mechanics with a production rate related directly to the volume. Compared with the classical theory of viscous flow [20,21], the production of inner energy due to dissipation is twice the value occurring in equation (4.13). Since, however, by (4.12) an additional production of inner energy at the inner boundaries is revealed, giving a contribution of the same amount as (4.13), we reason that the total production of inner energy is in accordance with classical theory.
The question arises as to whether the occurrence of discontinuous interfaces inside the fluid flow is an artefact of the model or if such phenomena really exist on a microscopic scale. Although a final answer to this question cannot be provided since knowledge about the processes occurring in a fluid flow on the micro-scale is restricted at the present time, it can be conjectured what kind of effect slight changes of the model may cause. At least the model established in this paper accurately recovers the physics on a macroscopic scale, and there are dissipative phenomena with entropy production at discontinuous surfaces known for a long time: we refer in particular to the classical theory of shock waves  where entropy production at discontinuous surfaces is provoked by a rapid compression of gas. Here, the discontinuous phenomena are not compression waves but can be considered as ‘slip waves’ which becomes apparent by (4.9).
Within this context, it is also of particular interest how the physics would be affected by a change of branch cut for the complex logarithm: in §(b) the standard branch cut along the negative real axis has been used, leading to values of between −iπ and iπ. One consequence of this is that over time the fluctuating forces (4.3) occurring in the equations of motion statistically result in zero by averaging. If an alternative branch cut for the complex logarithm is considered, say e.g. a cut at arguments of the complex number if a−π, the positions of the discontinuous surfaces are shifted and the values of the complex logarithm go from ia−iπ to ia+iπ. The first effect, the shift of the discontinuous surfaces, could be compensated for by applying the gauge transformation on the field of thermal excitation, whereas the second effect causes a change of the fluctuating forces (4.3) according to: The extra term on the right hand side vanishes in the limit ; however, for finite values of ω0 it gives an additional contribution to the equations of motion. By averaging again, this additional term does not tend to zero which is a valid physical argument for the use of the standard branch cut along the negative real axis (i.e. a=0).
Through the detailed analysis in §4, it has been demonstrated that the dynamics resulting from the Lagrangian (2.31) can self-consistently be interpreted as an extension of the classical theory towards processes beyond thermodynamic local equilibrium. A reproduction of the classical theory is reached by applying the limit for the relaxation rate to the equations of motion, a procedure that was not possible for the earlier suggested Lagrangian (2.4). In the following, some further indications are given in order to confirm the non-equilibrium assumption.
As already stated in §(a) in the context of the Lagrangian (2.4), a striking non-classical feature is the difference between the momentum density and the mass flux density, namely the quasi-momentum density p*=p−ϱu. In the context of the Lagrangian (2.31), the mass flux density is given via (A.18), whereas the momentum density is obtained as a canonical Noether observable , giving: . Hence, there is again a non-vanishing quasi-momentum density, 5.1 which in contrast to the quasi-momentum density (2.8) resulting from the Lagrangian (2.4) tends to zero for the limiting case . This is again in accordance with classical continuum theory.
Following the suggestion in previous work [6,12,13] that the quasi-momentum takes into account contributions to the momentum beyond the scope of the continuum hypothesis on a molecular scale, e.g. Brownian motion, a possible physical interpretation of the quasi-momentum density (2.8) is based on the elementary mechanism of viscosity, namely the transport of momentum by Brownian motion crosswise to the flow direction, see figure 6. The viscosity of a fluid is usually explained on a molecular scale by an exchange of particles between neighbouring fluid layers by Brownian motion of the molecules, by which a diffusion of momentum is induced. From the continuum viewpoint, the migrating particles responsible for the diffusive momentum flux are ‘quasi-particles’, associated with an additional contribution to the momentum density. Hence, the quasi-momentum density can also be considered as a measure for the deviation from thermodynamic equilibrium.
Note that the interpretation of the additional terms in the equations of motion as physical non-equilibrium effects on a microscopic scale is also consistent with the weak violation of the continuum hypothesis imposed by the discontinuous Lagrangian (2.31). Regardless, the violation of the continuum hypothesis vanishes in the limit for the relaxation rate: for increasing ω0, the discontinuities are decreasing and for very large ω0 they are physically reduced to fluctuations on a micro-scale in accordance with classical theory.
6. Conclusion and outlook
Based on an analogy between quantum mechanics and fluid mechanics, the formerly proposed Lagrangian (2.4) has been refined, leading to the discontinuous Lagrangian (2.31). By a careful analysis, it is proved that the dynamics resulting from Hamilton's principle based on (2.31) can consistently be interpreted as a generalization of the theory of viscous flow towards thermodynamic non-equilibrium, with a recovery of the classical Navier–Stokes equations and the balance of inner energy when applying the limit to the resulting equations of motion. As a striking feature, the application of the limit directly to the Lagrangian (2.31) fails. Hence, at least an indication is given that a variational formulation of viscous flow cannot be achieved using a continuous Lagrangian.
Although for large ω0 the discontinuities are physically reduced to fluctuations on a micro-scale, it would be of great interest for future work to explore further the dynamics of the discontinuous interfaces and the induced physical effects beyond thermodynamic equilibrium for various flow geometries.
As already mentioned in the Introduction, a striking feature for variational formulations in continuum mechanics is the necessity for the use of potentials in general . Since the Euler–Lagrange equations resulting from the two Lagrangians (2.4) and (2.31), the latter one is derived in §A.3 of appendix A, can also be interpreted as a first integral of the equations of motion, the use of potential fields seems also to be inevitable for the construction of first integrals of the equations of motion in fluid mechanics, see for example the attempts of He [28,29] for finding a potential representation of the velocity for two-dimensional incompressible and inviscid flow. For viscous flow, a first integral approach has been established using a generalized form of the Clebsch transformation , and based on the use of potentials [31,32], methodically different from the Clebsch transformation but with some interesting parallels to the approach given in this paper. The latter will be analysed in forthcoming papers.
By evaluating the dynamics induced by the Lagrangian (2.31), it has been demonstrated how Lagrange formalism applies to physical problems with discontinuities. Independent of the particular problem of viscous flow, the general formalism specified in §3 may also be a valuable mathematical tool for embedding various discontinuous phenomena into Lagrange formalism, like for example phase boundaries between non-mixable fluids, propagating shock waves in gaseous media, flame fronts, detonation shocks and also interfaces in solids like micro-cracks  and grain boundaries. Discontinuities also occur in some optimum control problems . The extended formalism can be used for finite-element simulations of such phenomena without the imperative of considering the related matching conditions explicitly, since they result automatically from the respective Lagrangian. It is therefore realistic to expect an improvement of, for example, numerical algorithms.
This is a pure theoretical work without experimental or numerical data. The analytical calculations performed are transparently given in the accompanying appendix.
F.M. supervised the bachelor thesis  revealing the failure of the conventional approach regarding the Poiseuille flow in a curved channel, worked out the parallels between the first integral approach given in earlier work [31,32] and the non-conventional approach given here, did the major part of the literature research and helped draft the manuscript; M.S. worked out the remaining benchmark tests, suggested the discontinuous Lagrangian (2.31), developed the extended formalism for discontinuous Lagrangians and drafted the manuscript. Both authors gave final approval for publication.
We declare we have no competing interests.
F.M. acknowledges the financial support from the Thomas Gessmann-Stiftung for his doctoral project. M.S. and F.M. acknowledge the support from the DFG (Deutsche Forschungsgemeinschaft), SCHO 767/6-1.
We thank N. Aksel, M. Oberlack, A. F. Cheviakov and A. Kluwick for helpful discussions.
Appendix A. Calculations
A.1 Symmetries and associated Noether balances of the Lagrangian (2.4)
In Scholle , a general analysis is provided concerning the analytical structure of Lagrangians in continuum theories fulfilling invariance with respect to the full Galilei group. In the same paper, a general scheme for Lagrangians is constructed. Using Noether's theorem, canonical formulae give rise to the identification of the relevant observable fields like mass density and flux density, momentum density, stress tensor, energy density and Poynting vector.
Subsequently, the analysis given in Scholle  is rigorously applied to the Lagrangian (2.4): for simultaneous invariance with respect to time and space translations and Galilei boosts, a collective symmetry criterion, the duality criterion, A 1 has to be fulfilled where are the conventional time derivative, the dual time derivative and the generating field, respectively, and A 2 and A 3 are the dual transformation and the corresponding infinitesimal generator. In the case of the present Lagrangian (2.4), the dual transformation takes the form A 4 fulfilling criterion (A.1) as required. One consequence of it is that the mass balance, A 5 is automatically fulfilled and ϱu is identified as the mass flux density. In the following, we skip the computation of canonical stress tensor, energy density and Poynting vector form Noether's theorem and focus on the canonical momentum density which takes the form  A 6 and need not be identical to the mass flux density ϱu. Since the dual transformation formula (A.4) contains a real ∇ζ-dependence, the relation A 7 is given with the quasi-momentum density: A 8 Hence, a striking feature of the Lagrangian (2.4) is a difference between mass flux density and momentum density which is in contrast to classical continuum mechanics.
A.2 Steady solution of (2.9)–(2.11) for Poiseuille flow
Considering the flow geometry u=u(y)ex and assuming time-independent thermasy ϑ, the evolution equation (2.11) takes the form A 9 which obviously implies the general solution (2.15). Considering (A.9), equations (2.9) read Taking the component in the y-direction, it follows that for the pressure (η=ϱ0ν) with integration function f2(x). By inserting this into the equation for the x-direction, we obtain Since the left hand side of the above equation depends on x only and the right hand on y only, both sides are constant, i.e. and with K>0. Hence, the pressure is evaluated as A 10 whereas the velocity profile has to fulfil the nonlinear ODE A 11 Via the substitution g(y)=1/u(y), equation (A.11) simplifies to the second-order ODE g′′−Kg2=0. After multiplication with g′ the latter becomes integrable, implying the first-order ODE g′2−2Kg3/3=C. The integration constant, C, is evaluated at the middle of the channel, y=h/2, where the velocity takes its maximum, u(h/2)=umax, and therefore u′(h/2)=0, implying C=−2K/(3u3 max). Hence, after re-substituting u=1/g, results as a nonlinear first-order ODE for the velocity profile, with the positive sign valid for 0≤y<h/2 and the negative sign for h/2<y≤h. Using separation of variables, the above ODE can be solved, with the solution for y≤h/2 implicitly given as A 12 and A 13 with Gaussian hypergeometric function 2F1.
A.3 Euler–Lagrange equations of the Lagrangian (2.31)
The Lagrangian (2.31) is based on the following fields: the velocity u, the three Clebsch variables Φ, α, β and the complex field of thermal excitation χ. Subsequently, the associated Euler–Lagrange equations are computed according to (3.3). First, variation with respect to the Clebsch variable Φ delivers the continuity equation A 14 As a consequence, the identity ∇⋅(ξu)=u⋅∇ξ, which is fulfilled for any field ξ, is considered subsequently. Next, variation with respect to the two remaining Clebsch variables α and β lead to the transport equations A 15 and A 16 after simple mathematical manipulation. By variation with respect to the velocity u, A 17 is obtained. Finally, the Euler–Lagrange equation related to the variation with respect to leads to the evolution equation A 18 for the thermal excitation. Variation with respect to delivers the complex conjugate of (A.19).
A.4 Equations of motion
The Euler–Lagrange equations (A.16)–(A.19) are a first integral of the equations of motion, i.e. the latter can be obtained from their derivations as follows: considering the identities and the Euler–Lagrange equations (A.15)–(A.19), the material time derivative of (A.18) reads A 19 where the pressure p is according to given as the Lagrangian evaluated for real processes, as is well known from the classical literature . Finally, a balance for the inner energy, is derived from (A.19) according to: A 20
A.5 Derivation of the production condition (3.9) using distributions
Within the theory of distributions, a very elegant way is provided by which the generalized formalism for discontinuous Lagrangians can be understood in terms of conventional Lagrange formalism. Under the same assumptions made in §3, by A 21 a continuous reference Lagrangian ℓc is defined, where is the Heaviside function giving 0 for x<0 and 1 for x≥0. Equation (A.22) defines a decomposition of the entire Lagrangian into a continuous part and a sum of discontinuities; the derivative of it with respect to φ gives A 22 where δ(x) denotes Dirac's delta function. Using this, the Euler–Lagrange expression ELN defined according to (3.3) with respect to ψN=φ can be defined across the entire domain, leading to As a consequence of the above decomposition, the first term in equation (3.5) for ψN=φ results in A 23 In order to evaluate the integrals at the right hand with delta functions, a representation of dV in terms of local coordinates, dV =dξ dS is used where dS is the surface element of the interface Sn and dξ is related to the direction perpendicular to Sn. By the substitution dφ=|∇φ| dξ, the identity
is obtained, giving (A.24) the more convenient form A 24 Finally, considering the identity (3.6) and the vanishing of δφ at t=t1,2, the variation of the action integral, (3.5), turns out to be A 25 which apparently implies, next to the well-known Euler–Lagrange equation with respect to φ, the production condition (3.9).
- Received June 29, 2016.
- Accepted January 4, 2017.
- © 2017 The Authors.
Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.