Should non-relativistic Navier Stokes Equations be modified so that they become pseudo-Lorentz invariant? Choking mass flow seems to reflect the fact that fluid momentum density has a maximum value (in stationary conditions) equal to $\rho_*  c_*$ where $\rho_*$ is the critical mass density and $c_*$ is the critical velocity, which is closely related to the speed of sound (see Landau-Lifchitz,"Fluid Mechanics",  section 83).
If this result also held in transient flows, would it not imply that the Navier-Stokes Equations should be modified so that they explicitly exhibited pseudo-Lorentzian symmetry (in momentum density) instead of Galilean symmetry? 
The momentum density field would then become explicitly causal.
The velocity field does not seem to make much sense if detached from the mass density field.
What I have in mind was the fact that, in a stationary fluid flow, from thermodynamic considerations, it can be proven that:
$$\frac{d(\rho v)}{dv}=\rho\big[1-\frac{v^2}{c^2}\big]$$
Landau and Lifhchitz, "Fluid Mechanics", section 83.
where $\rho$ is the mass density, $v$ the local velocity and $c$ the local speed of sound. 
This indicates that the momentum density does have a maximum, at least in stationary flows (where the local velocity is equal to the local speed of sound).
Since the momentum Navier-Stokes equations are Galilean invariant, a sufficiently large pressure gradient should make possible to attain larger momentum densities and this contradicts the previous equation. 
I was wondering if one should not enforce the maximum local momentum density as a pseudo-Lorentzian symmetry. (In incompressible flows, at least, I think that one does need to enforce this symmetry).
 A: There is an extended literature on relativistic dissipative hydrodynamics, modifying the Navier-Stokes equations to accommodate relativity. There are a number of inequivalent ways of doing this, only some of them leading to a causal theory. A paper by Van and Viro https://arxiv.org/abs/1109.0985 (one of many) gives lots of recent references. See also the book ''Rational extended thermodynamics'' by Müller and Ruggeri.
A: First, one should note that the expected Lorentz symmetry should be local rather than global. This could be understood by considering two different regions of the flow with identical physical quantities except the direction of the flow. Lorentz symmetry must be different for those points in a way that no global redifinition of variables could achieve. So we must be dealing with a point (and time) dependent, curved Lorentzian geometry, in other words OP's goals would be achieved by reformulating Navier–Stokes of Euler equations in terms of analog model of gravity.
Another point to note, is that full nonlinear dynamics of fluid (at least for arbitrary equations of state) likely would not allow formulation in terms of covariant (under new Lorentzian symmetry) objects. So as a reasonable minimum one could hope for 
a Lorentzian geometry to appear in linearized equations for the perturbations, although it is possible that a fuller nonlinear description would be present under certain conditions (equation of state, geometry of the flow).
And indeed something like this has already been done. In 1980 paper of Unruh a resemblance between black hole and sonic horizon of transsonic flow has been noted. Since then there has been a lot of work done on acoustic analog gravity. For a simple introduction have a look at papers [1, 2], and for extended survey of literature (and a lot of other analog models of gravity) a review [3].
From paper [2]:

If a fluid is barotropic and inviscid, and the flow is irrotational
  (though possibly time dependent) then the equation of motion for the
  velocity potential describing an acoustic disturbance is identical to
  the d'Alembertian equation of motion for a minimally coupled massless
  scalar field propagating in a $(3+1)$–dimensional Lorentzian
  geometry
  \begin{equation}
\Delta \psi \equiv 
{1\over\sqrt{-g}} 
\partial_\mu 
\left( \sqrt{-g} \; g^{\mu\nu} \; \partial_\nu \psi \right) = 0.
\end{equation}

The effective Lorentzian (acoustic) metric $g_{\mu\nu}$ would be
\begin{equation}
g_{\mu\nu} \equiv 
{\rho_0 \over  c}
\left[ \matrix{-(c^2-v_0^2)&\vdots&-v_0^j\cr
               \cdots\cdots\cdots\cdots&\cdot&\cdots\cdots\cr
        -v_0^i&\vdots&\delta_{ij}\cr } 
\right].        
\end{equation}
Equivalently, it can be expressed as
\begin{equation}
ds^2 \equiv g_{\mu\nu} \; dx^\mu \; dx^\nu =
{\rho_0\over c} 
\left[
- c^2 dt^2 + (dx^i - v_0^i \; dt) \; \delta_{ij} \; (dx^j - v_0^j \; dt )
\right].
\end{equation}
Note, that the above equation for perturbation is in terms of scalar velocity potential rather than for momentum density, but it is possible to construct appropriate covariant tensor objects and write covariant equations for them. 
Another thing to note, is that inclusion of viscosity will break the new Lorentz symmetry of the model (see section 12 of [2]).


*

*Visser, M. (1993). Acoustic propagation in fluids: an unexpected example of Lorentzian geometry. arXiv:gr-qc/9311028.

*Visser, M. (1998). Acoustic black holes: horizons, ergospheres and Hawking radiation. Classical and Quantum Gravity, 15(6), 1767, doi:10.1088/0264-9381/15/6/024, arXiv:gr-qc/9712010.

*Barcelo, C., Liberati, S., & Visser, M. (2011). Analogue gravity. Living reviews in relativity, 14(1), 3, doi:10.12942/lrr-2011-3.
A: Of course, as the velocity of the fluid approaches the speed of light, non-relativistic fluid dynamics should be replaced with a Lorentz invariant theory. However, for ordinary fluids (that excludes relativistic plasmas) the speed of sound is much less than the speed of light and ordinary Galilean invariant fluid dynamics should be perfectly good. 
The equation you quote does involve a few assumptions, steady flow, adiabatic flow, and no discontinuities (shocks). As you try to increase the pressure gradient these assumptions do not necessarily hold. In particular, the velocity field may steepen until you get a shock.
As Landau points out later in the book (sect 97), the equation you discuss
$$
 \frac{d(\rho v)}{dv} =   \rho \left( 1-\frac{v^2}{c_s^2}\right)
$$
has important consequences for the design of rocket engines. The greatest mass flow rate occurs at the nozzle constriction where $v\sim c_s$, and then the flow is supersonic in the expanding part of the nozzle. So this is an empirically verified equation, not one you would like to modify. 
Finally note that in the regime $v\sim c_s$ you cannot assume the flow to be incompressible (the speed of sound is the adiabativ compressibility of the fluid).  
