Derivation of Compressible Navier-Stokes Equations
Theorem
$$ \mathbf{u} = \mathbf{u} \left( t ; \mathbf{x} \right) = \left( u_{1} \left( t ; \mathbf{x} \right) , u_{2} \left( t ; \mathbf{x} \right) , u_{3} \left( t ; \mathbf{x} \right) \right) $$ In particular, in three-dimensional space, let the time $t$ and spatial coordinates $\mathbf{x} = \left( x_{1} , x_{2} , x_{3} \right)$ have the velocity field represented by the velocity vector as above. Similarly, $p : \mathbb{R}^{3} \to \mathbb{R}$ denotes the pressure $p = p \left( \mathbf{x} \right)$ at each location. If $\mathbf{u}$ is the velocity of a Newtonian fluid, it follows the following governing equation. $$ {\frac{ \partial \mathbf{u} }{ \partial t }} + \left( \mathbf{u} \cdot \nabla \right) \mathbf{u} = - \nabla w + \nu \nabla^{2} \mathbf{u} + \left( {\frac{ 1 }{ 3 }} \nu + \xi \right) \nabla \left( \nabla \cdot \mathbf{u} \right) + \mathbf{g} $$ Here $\nabla \cdot$ is the divergence, $\nu = \mu / \rho$ is the ratio of the dynamic viscosity $\mu$ to the density $\rho$, $\nabla w = \nabla p / \rho$ denotes thermodynamic work, $\xi = \zeta / \rho$ is the bulk kinematic viscosity (the ratio of the bulk viscosity $\zeta$ to the density), and $\mathbf{g}$ is the gravitational acceleration.
Explanation
At first glance the formulae may look intimidating by their length, but if you have followed the derivations of the Euler equations and the incompressible Navier–Stokes equations, you will see that this is merely a simple generalization obtained by adding terms one by one. $$ {\frac{ \partial \mathbf{u} }{ \partial t }} + \left( \mathbf{u} \cdot \nabla \right) \mathbf{u} = - \nabla w + \mathbf{g} $$ The Euler equations, as above, consider only pressure and gravity. $$ {\frac{ \partial \mathbf{u} }{ \partial t }} + \left( \mathbf{u} \cdot \nabla \right) \mathbf{u} = - \nabla w + \nu \nabla^{2} \mathbf{u} + \mathbf{g} $$ The Navier–Stokes equations add the viscous term $\nu \nabla^{2} \mathbf{u}$. $$ {\frac{ \partial \mathbf{u} }{ \partial t }} + \left( \mathbf{u} \cdot \nabla \right) \mathbf{u} = - \nabla w + \nu \nabla^{2} \mathbf{u} + \left( {\frac{ 1 }{ 3 }} \nu + \xi \right) \nabla \left( \nabla \cdot \mathbf{u} \right) + \mathbf{g} $$ The compressible Navier–Stokes equations cannot assume the incompressibility condition $\nabla \cdot \mathbf{u} = 0$ (equivalent), so the third term on the right-hand side simply remains.
Because compressible fluids usually refer to gases, it is acceptable to think of the incompressible Navier–Stokes equations as describing liquids and the compressible Navier–Stokes equations as describing gases.
Derivation
Divergence theorem: For a three-dimensional vector function $\mathbf{F}$ the following holds. $$ \int_{\mathcal{V}} \nabla \cdot \mathbf{F} dV = \oint_{\mathcal{S}} \mathbf{F} \cdot d \mathbf{S} $$ Here $\nabla \cdot \mathbf{F}$ is the divergence, $\int_{\mathcal{V}}$ is the volume integral, and $\oint_{\mathcal{S}}$ is the closed surface integral.
As in the derivation of the Navier–Stokes equations, we will apply the divergence theorem to the Cauchy stress tensor. If the sudden appearance of $\nabla \cdot \sigma$ is unclear, I recommend consulting the derivation of the Navier–Stokes equations, where this is explained more gently.
$$ \begin{align*} \mathbf{f} =& \oint_{\partial V} \sigma \cdot d S \\ =& \int_{V} \nabla \cdot \sigma d V \\ \implies {\frac{ d \mathbf{f} }{ d V }} =& \nabla \cdot \sigma \end{align*} $$
By the divergence theorem, adding $d \mathbf{f} / d V = \nabla \cdot \sigma$ to the right-hand side of $\rho D \mathbf{u} / D t$ yields the governing equation we want.
Cauchy stress tensor and Lamé parameters: In particular, in three-dimensional space, let the velocity field at time $t$ and spatial coordinates $\mathbf{x} = \left( x_{1} , x_{2} , x_{3} \right)$ be represented by the velocity vector as above. Suppose this fluid is a Newtonian fluid with viscosity and compressibility. The isotropic Cauchy stress tensor $\sigma$ can be expressed in terms of the symmetric gradient $\varepsilon$ as follows. $$ \sigma = - p I + 2 \mu \varepsilon + \lambda \tr \left( \varepsilon \right) I $$ Here the symmetric gradient $\varepsilon$ is defined as follows. $$ \varepsilon (\mathbf{u}) = {\frac{ 1 }{ 2 }} \left( \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{T} \right) $$
When both viscosity and compressibility are taken into account, the Cauchy stress tensor consists of three terms including pressure. To express these for $\mathbf{u}$, let us compute starting from $\tr \left( \varepsilon \right)$. $$ \tr \left( \varepsilon \right) = {\frac{ \partial u_{1} }{ \partial x_{1} }} + {\frac{ \partial u_{2} }{ \partial x_{2} }} + {\frac{ \partial u_{3} }{ \partial x_{3} }} = \nabla \cdot \mathbf{u} $$ Since $\varepsilon$ can be substituted directly by its definition, the vector form of $\sigma \left( \mathbf{u} \right)$ for the bulk viscosity $\zeta = \lambda + {\frac{ 2 }{ 3 }} \mu$ is as follows. $$ \begin{align*} & \sigma \left( \mathbf{u} \right) \\ =& - p I + \mu \left( \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{T} \right) + \lambda \left( \nabla \cdot \mathbf{u} \right) I \\ =& - p I + \lambda \left( \nabla \cdot \mathbf{u} \right) I + {\frac{ 2 }{ 3 }} \mu \left( \nabla \cdot \mathbf{u} \right) I - {\frac{ 2 }{ 3 }} \mu \left( \nabla \cdot \mathbf{u} \right) I + \mu \left( \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{T} \right) \\ =& - p I + \zeta \left( \nabla \cdot \mathbf{u} \right) I + \mu \left( \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{T} - {\frac{ 2 }{ 3 }} \left( \nabla \cdot \mathbf{u} \right) I \right) \end{align*} $$
Divergence of $\nabla \mathbf{u}$: $$ \nabla \cdot \nabla \mathbf{u} = \nabla^{2} \mathbf{u} $$ Divergence of $\left( \nabla \mathbf{u} \right)^{T}$: $$ \nabla \cdot \left( \nabla \mathbf{u} \right)^{T} = \nabla \left( \nabla \cdot \mathbf{u} \right) $$
Now let us expand $\nabla \cdot \sigma$ directly. $$ \begin{align*} & \rho {\frac{ D \mathbf{u} }{ D t }} - \rho \mathbf{g} \\ =& \nabla \cdot \sigma \\ =& \nabla \cdot \left( - p I + \zeta \left( \nabla \cdot \mathbf{u} \right) I + \mu \left[ \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{T} - {\frac{ 2 }{ 3 }} \left( \nabla \cdot \mathbf{u} \right) I \right] \right) \\ =& - \nabla p + \zeta \nabla \left( \nabla \cdot \mathbf{u} \right) + \mu \nabla \cdot \left[ \nabla \mathbf{u} + \left( \nabla \mathbf{u} \right)^{T} - {\frac{ 2 }{ 3 }} \left( \nabla \cdot \mathbf{u} \right) I \right] \\ =& - \nabla p + \zeta \nabla \left( \nabla \cdot \mathbf{u} \right) + \mu \left[ \nabla^{2} \mathbf{u} + \nabla \left( \nabla \cdot \mathbf{u} \right) - {\frac{ 2 }{ 3 }} \nabla \left( \nabla \cdot \mathbf{u} \right) \right] \\ =& - \nabla p + \zeta \nabla \left( \nabla \cdot \mathbf{u} \right) + \mu \left[ \nabla^{2} \mathbf{u} + {\frac{ 1 }{ 3 }} \nabla \left( \nabla \cdot \mathbf{u} \right) \right] \\ =& - \nabla p + \mu \nabla^{2} \mathbf{u} + {\frac{ 1 }{ 3 }} \mu \nabla \left( \nabla \cdot \mathbf{u} \right) + \zeta \nabla \left( \nabla \cdot \mathbf{u} \right) \\ =& - \nabla p + \mu \nabla^{2} \mathbf{u} + \left( {\frac{ 1 }{ 3 }} \mu + \zeta \right) \nabla \left( \nabla \cdot \mathbf{u} \right) \end{align*} $$
Finally, moving $\rho \mathbf{g}$ to the other side and dividing both sides by $\rho$ yields the equation we sought.
■