Fluid Dynamics in Curvilinear Coordinates without Fictitious Forces

Use of curvilinear coordinates is sometimes indicated by the inherent geometry of a fluid dynamics problem, but this introduces fictitious forces into the momentum equations that spoil strict conservative form. If one is willing to work in three dimensions, these fictitious forces can be eliminated by solving for rectangular (Cartesian) momentum components on a curvilinear mesh. A thoroughly geometric approach to fluid dynamics on spacetime demonstrates this transparently, while also giving insight into a greater unity of the relativistic and nonrelativistic cases than is usually appreciated.


Introduction
A fluid, like any material continuum, is defined and governed by tensor field equations on spacetime, a four-dimensional differentiable manifold [1,2]. The fluid's existence and kinematics are defined by the vanishing 4-divergence of a matter flux vector N, expressing conservation of the basic 'stuff' underlying the fluid (typically mass, or better, baryon number). The dynamics of the fluid is governed by Newton's second law and the first law of thermodynamics as applied to individual fluid elements. Given the 4-velocity field w defining the worldlines of fiducial observers, and in the absence of external forces and heating/cooling, these laws governing momentum and energy evolution can be combined in terms of the conservation law ∇ · S = 0 (2) for the relative energy-momentum flux, a (1, 1) tensor S [3,4]. The label "relative energymomentum flux" denotes the fact that the energy content of S includes only internal energy and bulk motion relative to fiducial observers w, allowing the same Equation (2) on spacetime to apply, regardless of whether causality is governed by absolute time (the "non-relativistic" world of Galilei and Newton) or by light cones (the "relativistic" world of Einstein and Minkowski) [3,4]. Even though there is no spacetime metric in the Galilei/Newton case, the 4-divergence still exists thanks to the presence of a 4-volume form (Levi-Civita tensor) E, through which defines the divergence of a vector field A; the right-hand side is the exterior derivative of the interior product of A with E.
To address an initial value problem, foliate spacetime into spacelike slices. For present purposes let spacetime be flat, a 4-dimensional affine space, either Galilei/Newton spacetime G equipped with absolute time [5][6][7][8][9][10][11][12][13] or Minkowski spacetime M equipped with a Lorentz metric [14]. In these cases we can choose the fiducial observers to be inertial observers, with worldlines given by the straight coordinate curves of a global time coordinate t, and 4-velocities w = ∂ t given by the natural basis vectors ∂ t tangent to this congruence of straight worldlines. On G the coordinate t is the absolute time, and the spacelike slices are already given as its level surfaces. On M the coordinate t is the proper time of the fiducial (and here, inertial) observers; and by virtue of the Lorentz metric, choose to specify that the spacelike slices be orthogonal to w, and thus level surfaces of t in this case as well. On both G and M these spacelike slices (S t ) are three-dimensional affine hyperplanes equipped with a Euclidean 3-metric γ. And in both cases the 4-volume form reads in which is the 3-volume form on the slices (S t ), expressed in terms of rectangular coordinates (x¯ı) = (x, y, z) or curvilinear coordinates u i = (u, v, w), with γ being the determinant of the matrix of curvilinear metric components. (Throughout this paper, latin indices denote coordinates, and also components of tensors, on the (S t ). Indices with an overbar are associated with rectangular coordinates, while unadorned indices denote curvilinear coordinates. Writing the (x¯ı) as functions of the u i , the elements of the matrix of curvilinear metric components γ ij are read off the line element dx 2 + dy 2 + dz 2 = γ ab du a du b . Note our use throughout of the Einstein summation convention; for readability, repeated summation indices a, b, c, . . . are taken from near the beginning of the alphabet, while i, j, k, . . . denote free indices.) This foliation of spacetime-here, into affine hyperplanes of G and M-effects a 3 + 1 decomposition of Equations (1) and (2) in terms of matter, momentum, and energy densities N, Π, and E respectively, and corresponding spatial fluxes F N , F Π , and F E tangent to S t , measured by the fiducial observers: in which the 3-divergences D · F are on S t , and √ γ has been taken to be independent of t. Expressing these in terms of component fields, highlights that while the matter and energy densities N and E are scalar fields, and their fluxes F N and F E are vector fields, the momentum density Π is a linear form, and its flux F Π is a (1, 1) tensor. For computations, these relations must be reduced to partial differential equations by further spelling out the 3-divergences on S t , yielding different expressions for different coordinate choices. Consider three ways to obtain the resulting equations, taking note of what happens to the suggestive conservative form of the original geometric expressions in terms of divergences.
First, under present assumptions, because S t is Euclidean (thanks to our focus here on flat G and M) we could begin with rectangular coordinates (x¯ı), with respect to which the covariant derivative operator Dī is simply a partial derivative ∂/∂x¯ı: ∂Πī ∂t

∂E ∂t
This strict conservative form can be translated into numerical methods (such as finite-volume discretization) that naturally handle discontinuities and reproduce global conservation to numerical precision. But under the coordinate transformation additional terms arise from derivatives of the position-dependent Jacobian factors, and this apparently spoils the conservative form.
Rather than proceed further with brute force, consider a second approach and make use of the rules of tensor analysis. In a coordinate basis, the covariant derivative of a tensor adds to the partial derivative an additional (sum of) "connection" terms for each tensor index. The gradient of a vector field yields the (1, 1) tensor field with components in which the connection coefficients Γ i jk (not tensor components in themselves!) are given by the Christoffel symbols in terms of derivatives of the metric components. (The (0, 2) metric tensor γ being nondegenerate, the matrix γ ij has matrix inverse γ ij gathering the components of the (2, 0) inverse metric tensor ← → γ .) The divergence of a vector field is the contraction of Equation (16). It turns out that Γ a ia = ∂ ln √ γ /∂u i , so that the divergence of a vector field can be expressed in the conservative form This is happily amenable to a structured grid in curvilinear coordinates, with a conservative finite-volume discretization of the divergence corresponding to the familiar definition encountered in elementary vector calculus: in which V ↔ is a finite cell volume, and the cell face areas A i and flux components F i are evaluated on the outer (i →) and inner (← i) faces in each dimesnion i. However, the momentum flux is a (1, 1) tensor field, and its covariant gradient has an additional term: Upon contraction to form the divergence D a F a i , the first connection term combines with the partial derivative as in the vector field case. And when F ij = F i a γ aj is symmetric-which is true for the momentum flux in fluid dynamics-there is some simplification in the second connection term, but it cannot be combined with the partial derivative. (This simplification is a consequence of choosing to solve for the covariant momentum components Π i rather than the contravariant ones Π i .) The fluid equations take the form ∂N ∂t ∂E ∂t The matter and energy equations are in strict conservative form. And the momentum equations might also be said to be in conservative form in a looser sense: it is a balance equation, with a "divergence" on the left-hand side and without derivatives of the fluid variables in the source terms, so that finite volume discretization can still handle discontinuities. But the terms on the right-hand side of Equation (22) constitute "fictitious forces": strictly speaking, curvilinear coordinates are non-inertial, and result in terms analogous to those that result from use of rotating or otherwise accelerated reference frames. Thus global conservation of curvilinear momentum components does not generally hold, and global conservation of rectangular momentum components will not be obtained to numerical precision. (Note in passing that in spherical and cylindrical coordinates, the azimuthal component of momentum is the component of angular momentum along the azimuthal axis, the equation for which is in strict conservative form. In cylindrical coordinates the equation for the linear momentum along the azimuthal axis is also in strict conservative form.) Use of a curvilinear coordinate mesh is sometimes indicated by the inherent geometry of a problem, even when three-dimensional phenomena preclude reduction in dimension because of an absence of axial or spherical symmetry. In such a case, is there any way to avoid the fictitious forces on the right-hand side of Equation (22)?
The purpose of this paper is to point out that the answer to this question is "yes", and to show it transparently, indeed almost instantly, with a third approach. If one is willing or needs to work in three dimensions anyway, solving for rectangular momentum components on a curvilinear coordinate patch allows for exploitation of spherical or cylindrical geometry without the fictitious forces on the right-hand side of Equation (22). This intuitively plausible result (mentioned without explanation or derivation in Ref. [15]) could of course be derived within the context of the first two approaches discussed in the preceding paragraphs. But it becomes particularly obvious when we step back from brute force coordinate transformations, or rules for tensor analysis on components, and consider the differential geometry behind those approaches, treating tensors as geometric objects [14,16].

Mixed basis for the momentum flux
The point is that a tensor field is not just its component functions (even though these uniquely determine it). Components are merely expansion factors appearing when a tensor is expressed in terms of a basis-or, for fields, component functions appear when a tensor field is expressed in terms of a smoothly varying basis field. The covariant derivative makes it possible to compare tensors at neighboring points of a manifold by taking account not only of the variation of component functions, but also the variation of the tensor basis field elements in terms of which the tensor field is expanded.
In particular, a (1, 1) tensor field is given by in terms of a basis field e i ⊗ e j formed, via the tensor product, from some basis field (e i ) of vectors and some basis field e i of linear forms. The tensor field itself-a geometric object-is the same, regardless of whether one chooses (for example) the natural vector basis (∂ x¯ı ) and its dual linear form basis (dx¯ı) of rectangular coordinates, or the natural bases (∂ u i ) and du i of curvilinear coordinates: We may be used to treating coordinate transformations as an all-or-nothing affair, for instance a choice between the two expressions in Equation (25), to give equations in which all indices correspond to a particular choice of coordinates. In fact, however, nothing prevents us from using the mixed curvilinear/rectangular basis on a coordinate patch with curvilinear coordinates u i . The covariant derivative obeys the Leibniz rule for derivatives of products, so that Rules for tensor analysis notwithstanding, the component functions F i are actually scalar fields, for which the covariant derivative coincides with the partial derivative: Moreover the derivatives of the basis vector and linear form fields are precisely what give the connection coefficients in Equation (20): But rectangular coordinate basis fields on Euclidean space are constant; any variation vanishes, including in particular D j dx¯ı = 0 (30) appearing in the third term of Equation (27). Therefore The divergence-the contraction on the first and third slots-takes the form since Recall also that Γ a ia = ∂ ln √ γ /∂u i , as noted previously. Thus Equation (32) shows that the fictitious forces on the right-hand side of Equation (  22), which stem from the connection term on the covariant index of a (1, 1) tensor F, vanish when F is expanded in terms of the mixed basis of Equation (26). Then, on flat spacetimes G and M and in the absence of external forces and heating/cooling, all the fluid equations ∂N ∂t ∂Πī ∂t ∂E ∂t are in strict conservative form. Again, it is the rectangular momentum component fields Πī being solved for, from the expansion Π = Πā dxā, rather than the curvilinear component field from the alternative expansion Π = Π a du a .

Fluid dynamics equations in curvilinear and mixed bases
By way of concrete example, consider the fluid dynamics equations on M and G, specializing to a perfect fluid, in curvilinear coordinates general enough to encompass the rectangular, cylindrical, and spherical cases, before and after eliminating the fictitious forces associated with these coordinate choices. Before arriving at partial differential equations in particular coordinates it is necessary to give a 3 + 1 decomposition of the spacetime tensor Equations (1) and (2) satisfied by the matter and relative energy-momentum 4-fluxes, taking account also of the necessity to find expressions in terms of quantities measured by a comoving observer in order to apply constitutive equations (in this case, an equation of state). For further background, and on allowance for heat flow and viscosity, see Refs. [3,4] (beware however some notational changes).

3 + 1 decomposition of the matter and relative energy-momentum fluxes
Recall the foliation of affine spacetimes M and G introduced in Section 1. Fiducial (and here, inertial) observers have straight worldlines with 4-velocity vector field w = ∂ t threading spacelike slices (S t ), affine hyperplanes, the level surfaces of global time coordinate t. Introduce now the the linear form t = ∇t = dt associated with the (S t ). Because dt is dual while for a vector a tangent to S t , a · t = t · a = 0.
(Note that in this paper the dot, ·, only denotes contraction and never an inner product.) We have already introduced γ as the 3-metric on S t . It is consistent to use the same notation to denote the related projection tensor from M or G to S t . In particular we need the related projection tensor for a tangent to S t and the linear form the metric dual of a on S t . The conservative formulations we desire follow from decomposing vector fields on M or G into pieces parallel to w and tangent to S t . Under the geometry described above, a vector field A so decomposed, We so decompose the matter flux N and relative energy-momentum flux S. Let the flow of matter-taken to be baryon number-define a field of comoving observers. That is, where n is the baryon number density measured by the comoving observers, and U is their 4-velocity field. From the perspective of the fiducial observers, where v is tangent to S t and is the 3-velocity of the fluid measured by the fiducial observers. On M, the Lorentz factor Λ v = 1 − γ(v, v)/c 2 follows from the normalization g(U, U) = −c 2 given by the Lorentz metric g. Comparing with the decomposition we find and F N = Nv on M and G.
Therefore, in accord with Equation (44), ∇ · N = 0 becomes with however the differing baryon densities N measured by fiducial observers on M and G as given by Equation (48). The flow of relative energy-momentum calls for some explanation [4]. It is a (1, 1) tensor expressed in terms of two primary pieces. The first term in Equation (51) is the relative 4-momentum per baryon P, a linear form, carried by the baryon flux nU. To begin to understand P, consider first the existence on both M and G of an inertia flux vector I = mU, where m is the mass per baryon. The inertia flux has timelike component Λ v m on M and m on G, the inertia per baryon measured by a fiducial observer. The metric dual of I on M with respect to g, the linear form g · I, is the traditional relativistic 4-momentum per baryon, with timelike component −Λ v mc 2 , the (negative of the) mass-energy per baryon as measured by a fiducial observer. But this timelike component becomes meaningless as c → ∞: unlike the inertia flux, the traditional relativistic 4-momentum makes no sense on G, which allows no mass-energy equivalence. What does make sense on both M and G is to define a "relative" 4-momentum per baryon P, a linear form whose timelike component is the kinetic energy measured by the fiducial observers. On M it can be expressed as g · m(U − w), yielding with the definition on G taken to be the c → ∞ limit of the expression on M. Recall that v = γ · v is a linear form on S t , the metric dual of the 3-velocity v, which is tangent to S t . The second term in Equation (51), the (negative of the) 4-stress Σ, a (1, 1) tensor, arises from considering each fluid element to be an infinitesimal thermodynamic system in its own right (this tensor can also be arrived at using a variational approach [17]). For physical interpretation it is once again convenient to begin with a raised index, considering the (2, 0) tensor − → Σ = Σ · ← → g on M and the (negative of) inertia-momentum flux it represents. For a perfect fluid, The first term is the flux of inertia per baryon due to internal energy density e, measured by comoving observers (and with no heat flux out of fluid elements), carried by the baryon flux nU. In the second term, the projection tensor orthogonal to U enforces the spatial stress as seen by comoving observers to be isotropic and given by a pressure p. Lowering the second index gives with the definition on G following once again as the c → ∞ limit of the expression on M. Note that S must be a (1, 1) tensor in order combine energy and momentum on G (even though it can only be kinetic plus internal energy-that is, macroscopic plus microscopic kinetic energy) [4]. A contravariant index is needed to have a flow on spacetime, and a divergence of that flow. But a covariant index is also needed, in order for energy in a timelike component to survive the c → ∞ limit. This seems to be a reflection of absolute inertia corresponding to absolute time on G. In any case, that 4-momentum be consistently regarded as naturally a linear form, while 4-velocity is a vector, is consistent with the Hamiltonian perspective of (4-)momentum being conjugate to (4-)position.
Putting Equations (52) and (54) in Equation (51), the relative energy-momentum flux can be expressed as where and In a moment it will become clear that the linear form Π is the momentum density of the fluid measured by fiducial observers, and that the scalar field E is the relative energy density of the fluid measured by fiducial observers. Because S is a rank-2 tensor, its 3 + 1 decomposition must proceed on both indices. Begin with the second (covariant) index, using S as given in Equation (55). Isolate the momentum flux M, also a (1, 1) tensor field, as the spatial projection of S according to fiducial observers: where the terms in parentheses in the last line are the momentum 3-flux on S t . Isolate the relative energy flux E, a vector field, as the projection of S along w: For our inertial observers on flat spacetimes M and G, the covariant derivatives of w and ← − γ vanish. Therefore ∇ · S = 0 breaks up into ∇ · M = 0 and ∇ · E = 0; and according to Equation (44), these give with however differing momentum and relative energy densities Π and E measured by fiducial observers on M and G, as given by Equations (56) and (57).
For cylindrical coordinates ( , z, φ) given by (x, y, z) = ( cos φ, sin φ, z), the inverse Jacobian is The dependence of the flux components on θ and φ means that computations must be performed in three dimensions: reduction of dimension under conditions of axial or spherical symmetry is not an option when solving for the rectangular components of momentum.

Conclusion
Some fluid dynamics problems are couched most naturally in terms of spherical or cylindrical coordinates, even when three-dimensional phenomena are of interest in a generally spherical or cylindrical geometric scenario. When one needs or is willing to work in three dimensions anyway, the fictitious forces arising from curvilinear coordinates can be eliminated by solving for rectangular momentum components on a curvilinear mesh. This allows a strictly conservative form of the fluid dynamics equations to be restored on flat spacetime and in the absence of external forces and heating/cooling. However, angular dependence is thereby introduced into the momentum flux: the advected rectangular momentum components do not exhibit spherical or cylindrical symmetry, and the stress is transformed by a Jacobian matrix that explicitly introduces angular dependence. (Thus the necessity to work in three dimensions if this option is to be exercised.) But fictitious force terms that grow large near coordinate singularities are absorbed into the flux, where they are partially ameliorated when multiplied by the metric determinant appearing in the divergence.
When spherical-or cylindrical-style coordinates are used in the curved spacetime of general relativity, this approach can still eliminate the source terms associated with those coordinates that may cause trouble near coordinate singularities, even as more physically relevant source terms representing gravitation remain.
The possibility of eliminating fictitious forces arising from curvilinear coordinates is made transparent by an approach and notation that treat tensors as geometric objects. The thoroughly geometric approach modeled here will also improve the transparency of formulations involving curved spacetime or otherwise accelerated, i.e. Lagrange or arbitrary-Euler-Lagrange, observers. (The difference from the present focus on inertial fiducial observers in flat spacetime will be that the 4-velocity and spatial projection tensor fields associated with accelerated observers will not have vanishing covariant derivative.) This geometric approach also highlights a greater unity between relativistic and non-relativistic fluid dynamics than is usually appreciated.