A note on reverse engineering the Navier-Stokes equations

The Navier-Stokes equations are the fundamental equations of fluid mechanics, analogous to, say, Maxwell’s equations in the case of electromagnetism. They are important for applications in science and engineering but they are difficult to solve analytically so real-life applications rely almost exclusively on computer-aided methods. In fact, the equations are still not yet fully understood mathematically. Some basic questions about the existence and nature of possible solutions remain unanswered. Because of their importance, the Clay Mathematics Institute has offered a prize of one million dollars to anyone who can clarify some specific fundamental questions about them (see the CMI web page for details).

The equations are named after Claude-Louis Navier (1785-1836) and George Gabriel Stokes (1819-1903) who worked on their development independently. Some other leading contemporary mathematicians were also involved and I was intrigued to learn that all of these mathematicians used Euler’s equations of motion as a starting point in their derivations of the Navier-Stokes equations. Euler’s equations had been derived much earlier in 1757 by the great mathematician Leonhard Euler (1707-1783). In tensor notation, using the Einstein summation convention, the Euler equations take the form

$\rho g^i - \frac{\partial p}{\partial x^i} = \rho \big(\frac{\partial v^i}{\partial t} + v^j \frac{\partial v^i}{\partial x^j}\big)$

I explored a derivation of these equations from Newton’s second law in a previous post. To better understand how the two sets of equations are related, it occurred to me that it might be an interesting exercise, and probably not too difficult, to try to reverse engineer’ the Navier-Stokes equations using Euler’s equations as a starting point. I want to briefly record my exploration of this idea in the present note.

I am interested in the incompressible fluid version of the Navier-Stokes equations which can be written in tensor form, again using the Einstein summation convention, as

$\rho g^i - \frac{\partial p}{\partial x^i} + \mu \frac{\partial}{\partial x^j}\big(\frac{\partial v^i}{\partial x^j}\big) = \rho \big(\frac{\partial v^i}{\partial t} + v^j \frac{\partial v^i}{\partial x^j}\big)$

where $\mu > 0$ is a viscosity coefficient specific to the fluid in question. Putting the Euler equations and the Navier-Stokes equations side by side like this makes it clear that the key difference between them is the addition of the viscosity-related terms on the left-hand side of the Navier-Stokes equations. To clarify things a bit more it is helpful to see the equations laid out fully. The Euler equations constitute a three-equation system which looks like this:

$\rho g^1 - \frac{\partial p}{\partial x^1} = \rho \big(\frac{\partial v^1}{\partial t} + v^1 \frac{\partial v^1}{\partial x^1} + v^2 \frac{\partial v^1}{\partial x^2} + v^3 \frac{\partial v^1}{\partial x^3}\big)$

$\rho g^2 - \frac{\partial p}{\partial x^2} = \rho \big(\frac{\partial v^2}{\partial t} + v^1 \frac{\partial v^2}{\partial x^1} + v^2 \frac{\partial v^2}{\partial x^2} + v^3 \frac{\partial v^2}{\partial x^3}\big)$

$\rho g^3 - \frac{\partial p}{\partial x^3} = \rho \big(\frac{\partial v^3}{\partial t} + v^1 \frac{\partial v^3}{\partial x^1} + v^2 \frac{\partial v^3}{\partial x^2} + v^3 \frac{\partial v^3}{\partial x^3}\big)$

The Navier-Stokes equations also constitute a three-equation system which looks like this:

$\rho g^1 - \frac{\partial p}{\partial x^1} + \mu \bigg(\frac{\partial}{\partial x^1}\big(\frac{\partial v^1}{\partial x^1}\big) + \frac{\partial}{\partial x^2}\big(\frac{\partial v^1}{\partial x^2}\big) + \frac{\partial}{\partial x^3}\big(\frac{\partial v^1}{\partial x^3}\big)\bigg) = \rho \big(\frac{\partial v^1}{\partial t} + v^1 \frac{\partial v^1}{\partial x^1} + v^2 \frac{\partial v^1}{\partial x^2} + v^3 \frac{\partial v^1}{\partial x^3}\big)$

$\rho g^2 - \frac{\partial p}{\partial x^2} + \mu \bigg(\frac{\partial}{\partial x^1}\big(\frac{\partial v^2}{\partial x^1}\big) + \frac{\partial}{\partial x^2}\big(\frac{\partial v^2}{\partial x^2}\big) + \frac{\partial}{\partial x^3}\big(\frac{\partial v^2}{\partial x^3}\big)\bigg) = \rho \big(\frac{\partial v^2}{\partial t} + v^1 \frac{\partial v^2}{\partial x^1} + v^2 \frac{\partial v^2}{\partial x^2} + v^3 \frac{\partial v^2}{\partial x^3}\big)$

$\rho g^3 - \frac{\partial p}{\partial x^3} + \mu \bigg(\frac{\partial}{\partial x^1}\big(\frac{\partial v^3}{\partial x^1}\big) + \frac{\partial}{\partial x^2}\big(\frac{\partial v^3}{\partial x^2}\big) + \frac{\partial}{\partial x^3}\big(\frac{\partial v^3}{\partial x^3}\big)\bigg) = \rho \big(\frac{\partial v^3}{\partial t} + v^1 \frac{\partial v^3}{\partial x^1} + v^2 \frac{\partial v^3}{\partial x^2} + v^3 \frac{\partial v^3}{\partial x^3}\big)$

Looking at the pattern of the indices appearing in the viscosity-related terms on the left-hand side of the Navier-Stokes equations, one is immediately reminded of the indices identifying the positions of the elements of a matrix and we can in fact arrange the viscosity-related terms in matrix form as

$\begin{bmatrix} \frac{\partial}{\partial x^1}\big(\frac{\partial v^1}{\partial x^1}\big) & \frac{\partial}{\partial x^2}\big(\frac{\partial v^1}{\partial x^2}\big) & \frac{\partial}{\partial x^3}\big(\frac{\partial v^1}{\partial x^3}\big)\\ \ \\ \frac{\partial}{\partial x^1}\big(\frac{\partial v^2}{\partial x^1}\big) & \frac{\partial}{\partial x^2}\big(\frac{\partial v^2}{\partial x^2}\big) & \frac{\partial}{\partial x^3}\big(\frac{\partial v^2}{\partial x^3}\big)\\ \ \\ \frac{\partial}{\partial x^1}\big(\frac{\partial v^3}{\partial x^1}\big) & \frac{\partial}{\partial x^2}\big(\frac{\partial v^3}{\partial x^2}\big) & \frac{\partial}{\partial x^3}\big(\frac{\partial v^3}{\partial x^3}\big)\end{bmatrix}$

From the point of view of reverse engineering the Navier-Stokes equations this is highly suggestive because in deriving the Euler equations we encountered one key matrix, the rank-2 stress tensor $\tau^{ij}$. Could the above matrix be related to the $3 \times 3$ stress tensor? The answer is yes, and in fact we have that

$\begin{bmatrix} \frac{\partial}{\partial x^1}\big(\frac{\partial v^1}{\partial x^1}\big) & \frac{\partial}{\partial x^2}\big(\frac{\partial v^1}{\partial x^2}\big) & \frac{\partial}{\partial x^3}\big(\frac{\partial v^1}{\partial x^3}\big)\\ \ \\ \frac{\partial}{\partial x^1}\big(\frac{\partial v^2}{\partial x^1}\big) & \frac{\partial}{\partial x^2}\big(\frac{\partial v^2}{\partial x^2}\big) & \frac{\partial}{\partial x^3}\big(\frac{\partial v^2}{\partial x^3}\big)\\ \ \\ \frac{\partial}{\partial x^1}\big(\frac{\partial v^3}{\partial x^1}\big) & \frac{\partial}{\partial x^2}\big(\frac{\partial v^3}{\partial x^2}\big) & \frac{\partial}{\partial x^3}\big(\frac{\partial v^3}{\partial x^3}\big)\end{bmatrix} = \begin{bmatrix} \frac{\partial \tau^{11}}{\partial x^1} & \frac{\partial \tau^{12}}{\partial x^2} & \frac{\partial \tau^{13}}{\partial x^3}\\ \ \\ \frac{\partial \tau^{21}}{\partial x^1} & \frac{\partial \tau^{22}}{\partial x^2} & \frac{\partial \tau^{23}}{\partial x^3}\\ \ \\ \frac{\partial \tau^{31}}{\partial x^1} & \frac{\partial \tau^{32}}{\partial x^2} & \frac{\partial \tau^{33}}{\partial x^3} \end{bmatrix}$

In deriving the Euler equations what we did was to assume that the fluid was friction-free so that there were no shear forces acting on any of the faces of a given fluid volume element. This means that the fluid was being modelled as being inviscid, i.e., having no viscosity. As a result, we kept only the diagonal elements of the stress tensor and re-interpreted these as pressures. Approximating the net pressure in each direction on a fluid volume element using a Taylor series expansion we were then led to include only the first-order partials of the diagonal terms of the stress tensor (interpreted as the first-order partials of pressure) in the Euler equations.

In the case of the Navier-Stokes equations we are no longer assuming that the fluid is inviscid so we are including all the first-order partials of the stress tensor in the equations. Thus, the Navier-Stokes equations could be written as

$\rho g^1 - \frac{\partial p}{\partial x^1} + \mu \big(\frac{\partial \tau^{11}}{\partial x^1} + \frac{\partial \tau^{12}}{\partial x^2} + \frac{\partial \tau^{13}}{\partial x^3}\big) = \rho \big(\frac{\partial v^1}{\partial t} + v^1 \frac{\partial v^1}{\partial x^1} + v^2 \frac{\partial v^1}{\partial x^2} + v^3 \frac{\partial v^1}{\partial x^3}\big)$

$\rho g^2 - \frac{\partial p}{\partial x^2} + \mu \big(\frac{\partial \tau^{21}}{\partial x^1} + \frac{\partial \tau^{22}}{\partial x^2} + \frac{\partial \tau^{23}}{\partial x^3}\big) = \rho \big(\frac{\partial v^2}{\partial t} + v^1 \frac{\partial v^2}{\partial x^1} + v^2 \frac{\partial v^2}{\partial x^2} + v^3 \frac{\partial v^2}{\partial x^3}\big)$

$\rho g^3 - \frac{\partial p}{\partial x^3} + \mu \big(\frac{\partial \tau^{31}}{\partial x^1} + \frac{\partial \tau^{32}}{\partial x^2} + \frac{\partial \tau^{33}}{\partial x^3}\big) = \rho \big(\frac{\partial v^3}{\partial t} + v^1 \frac{\partial v^3}{\partial x^1} + v^2 \frac{\partial v^3}{\partial x^2} + v^3 \frac{\partial v^3}{\partial x^3}\big)$

Detailed consideration of the forces acting on an infinitesimal volume element in the case of a viscous fluid lead to expressions for the stress tensor components as functions of the first-order partials of the components of the velocity vector of the form

$\tau^{ij} = \tau^{ji} = \mu \big(\frac{\partial v^i}{\partial x^j} + \frac{\partial v^j}{\partial x^i}\big)$

The first-order partials of these are then of the form

$\frac{\partial \tau^{ij}}{\partial x^j} = \mu \frac{\partial }{\partial x^j}\big(\frac{\partial v^i}{\partial x^j}\big) + \mu \frac{\partial }{\partial x^j}\big(\frac{\partial v^j}{\partial x^i}\big)$

But note that with incompressible flow the mass density of the fluid is a constant so the divergence of the velocity vanishes (this is implied directly by the continuity equation with constant mass density). Therefore the second term on the right-hand side is zero since

$\frac{\partial }{\partial x^j}\big(\frac{\partial v^j}{\partial x^i}\big) = \frac{\partial }{\partial x^i}\big(\frac{\partial v^j}{\partial x^j}\big)$

$= \frac{\partial }{\partial x^i}\big(\frac{\partial v^1}{\partial x^1}\big) + \frac{\partial }{\partial x^i}\big(\frac{\partial v^2}{\partial x^2}\big) + \frac{\partial }{\partial x^i}\big(\frac{\partial v^3}{\partial x^3}\big) = \frac{\partial }{\partial x^i} \text{div}(\vec{v}) = 0$

Thus we are left with

$\frac{\partial \tau^{ij}}{\partial x^j} = \mu \frac{\partial }{\partial x^j}\big(\frac{\partial v^i}{\partial x^j}\big)$

which are the extra terms appearing on the left-hand side of the Navier-Stokes equations.

A fluid-mechanical visualisation of the quantum-mechanical continuity equation

The concept of a probability current is useful in quantum mechanics for analysing quantum scattering and tunnelling phenomena, among other things. However, I have noticed that the same rather abstract and non-visual approach to introducing probability currents is repeated almost verbatim in every textbook (also see, e.g., this Wikipedia article). The standard approach essentially involves defining a probability current from the outset as

$\vec{j} = \frac{i \hbar}{2m}(\Psi \nabla \Psi^{*} - \Psi^{*} \nabla \Psi)$

and then using Schrödinger’s equation to show that this satisfies a fluid-like continuity equation of the form

$\frac{\partial \rho}{\partial t} + \nabla \cdot \vec{j} = 0$

with

$\rho \equiv \Psi^{*} \Psi$

In the present note I want to briefly explore a more intuitive and visual approach involving a model of the actual flow of a probability fluid’. I want to begin with a fluid-mechanical model and then obtain the standard expression for the quantum-mechanical continuity equation from this, rather than starting out with an abstract definition of the probability current and then showing that this satisfies a continuity equation. The essential problem one faces when trying to do this is that although in classical mechanics the position $\vec{r}(t)$ of a point particle and its velocity $\vec{v}(t) = d\vec{r}(t)/dt$ are well defined, this is not the case in conventional quantum mechanics. Quantum mechanics is done probabilistically, treating a particle as a wave packet such that the square of the amplitude of the corresponding wave function acts as a probability density which can be used to measure the probability that the particle will occupy a particular region of space at a particular time. It is not possible to say definitively where a particular particle will be at a particular time in quantum mechanics, which makes it difficult to apply the conventional deterministic equations of fluid mechanics.

A huge specialist literature on quantum hydrodynamics has in fact arisen which tries to circumvent this problem in a number of ways. A standard reference is Wyatt, R. E., 2005, Quantum Dynamics with Trajectories: Introduction to Quantum Hydrodynamics (Springer). The route that a large part of this literature has taken is intriguing because it is based on Bohmian mechanics, an approach to quantum mechanics developed by David Bohm in 1952 which is regarded by most mainstream physicists today as unconventional. The key feature of the Bohmian mechanics approach is that classical-like particle trajectories are possible. Using this approach one can obtain Newtonian-like equations of motion analogous to those in conventional fluid mechanics and this is how this particular literature seems to have chosen to treat quantum particle trajectories in a fluid-like way. Attempts have also been made to introduce mathematically equivalent approaches, but defined within conventional quantum mechanics (see, e.g., Brandt, S. et al, 1998, Quantile motion and tunneling, Physics Letters A, Volume 249, Issue 4, pp. 265-270).

In the present note I am not looking to solve any elaborate problems so I will simply consider a free quantum wave packet which is not acted upon by any forces and try to visualise probability currents and the quantum continuity equation in a fluid-like way by using the momentum vector operator $\hat{\vec{p}}$ to characterise the velocity of the particle. I will then show that the probability current obtained in this fluid-mechanical model is the same as the one defined abstractly in textbooks.

In quantum mechanics, calculations are done using operators to represent observables. Every possible observable that one might be interested in for the purposes of experiment has a corresponding operator which the mathematics of quantum mechanics can work on to produce predictions. The key operator for the purposes of the present note is the momentum vector operator

$\hat{\vec{p}} = -i \hbar \nabla$

which is the quantum mechanical analogue of the classical momentum vector

$\vec{p} = m \vec{v}$

The key idea for the present note is to regard the velocity vector of the quantum particle as being represented by the operator

$\frac{\hat{\vec{p}}}{m}$

by analogy with the classical velocity vector which can be obtained as

$\vec{v} = \frac{\vec{p}}{m}$

We will imagine the total probability mass

$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} \Psi^{*} \Psi \text{d}V = 1$

as a fluid in steady flow throughout the whole of space and obeying mass conservation. The fact that the flow is steady reflects the fact that there are no forces acting on the quantum particle in this model, so we must have

$\frac{\partial }{\partial t}\big[\frac{\hat{\vec{p}}}{m}\big] = 0$

The velocity can vary from point to point in the probability fluid but at any given point it cannot be varying over time.

In a classical fluid we have a mass density per unit volume $\rho$ and we regard the velocity vector $\vec{v}$ as a volumetric flow rate per unit area, i.e., the volume of fluid that would pass through a unit area per unit time. Then $\rho \vec{v}$ is the mass flow rate per unit area, i.e., the mass of fluid that would pass through a unit area per unit time. In quantum mechanics we can regard the probability mass density per unit volume $\rho \equiv \Psi^{*} \Psi$ as analogous to the mass density of a classical fluid. We can interpret $\frac{\hat{\vec{p}}}{m}$ as the volumetric flow rate per unit area, i.e., the volume of probability fluid that would pass through a unit area per unit time. When doing probability calculations with quantum mechanical operators we usually sandwich’ the operator between $\Psi^{*}$ and $\Psi$, so following that approach here we can define the probability current density as

$\Psi^{*} \frac{\hat{\vec{p}}}{m} \Psi$

This is to be interpreted as the probability mass flow rate per unit area, i.e., the amount of probability mass that would pass through a unit area per unit time, analogous to $\vec{j} = \rho \vec{v}$ in the classical case. To see how close the analogy is, suppose the quantum wave function is that of a plane wave

$\Psi(x, y, z, t) = A\mathrm{e}^{i(\vec{k} \cdot \vec{r} - \omega t)}$

Then

$\Psi^{*} \frac{\hat{\vec{p}}}{m} \Psi = A \mathrm{e}^{-i(\vec{k} \cdot \vec{r} - \omega t)} \frac{(-i \hbar)}{m} \nabla A \mathrm{e}^{i(\vec{k} \cdot \vec{r} - \omega t)}$

$= A \mathrm{e}^{-i(\vec{k} \cdot \vec{r} - \omega t)} \frac{(-i \hbar)}{m} A \ i \ \vec{k} \ \mathrm{e}^{i(\vec{k} \cdot \vec{r} - \omega t)}$

$= A^2 \frac{\hbar \vec{k}}{m}$

$= \rho \vec{v}$

which looks just like the mass flow rate in the classical case with $\rho = \Psi^{*} \Psi$ and $\vec{v} \equiv \frac{\hbar \vec{k}}{m}$. Note that in this example the probability current density formula we are using, namely $\Psi^{*} \frac{\hat{\vec{p}}}{m} \Psi$, turned out to be real-valued. Unfortunately this will not always be the case. Since the probability current vector must always be real-valued, the fluid-mechanical model in the present note will only be applicable in cases when this is true for the formula $\Psi^{*} \frac{\hat{\vec{p}}}{m} \Psi$.

As in classical fluid mechanics, a continuity equation can now be derived by considering the net outflow of probability mass from an infinitesimal fluid element of volume $\mathrm{d}V \equiv \mathrm{d} x \mathrm{d} y \mathrm{d} z$.

Considering only the $y$-component for the moment, we see from the diagram that on the left-hand side we have the probability mass flow rate coming into the volume element through the left-hand face. The mass flow rate coming out of the fluid element through the right-hand face can be approximated using a Taylor series expansion as being equal to the mass flow rate through the left-hand face plus a differential adjustment based on the gradient of the probability current density and the length $\mathrm{d} y$. The net probability mass flow rate in the $y$-direction is then obtained by subtracting the left-hand term from the right-hand term to get

$\frac{\partial }{\partial y} \big(\Psi^{*}\frac{\hat{\vec{p}}}{m} \Psi\big) \mathrm{d}V$

Using similar arguments for the $x$ and $z$-directions, the net mass flow rate out of the fluid element in all three directions is then

$\nabla \cdot \big(\Psi^{*}\frac{\hat{\vec{p}}}{m} \Psi \big) \mathrm{d} V$

Now, the probability mass inside the fluid element is $\rho \mathrm{d} V$ where $\rho = \Psi^{*} \Psi$ and if there is a net outflow of probability fluid this mass will be decreasing at the rate

$- \frac{\partial \rho}{\partial t} \mathrm{d} V$

Equating the two expressions and dividing through by the volume of the fluid element we get the equation of continuity

$\frac{\partial \rho}{\partial t} + \nabla \cdot \big(\Psi^{*}\frac{\hat{\vec{p}}}{m} \Psi \big) = 0$

What I want to do now is show that if we work out $\Psi^{*}\frac{\hat{\vec{p}}}{m} \Psi$ we will get the same formula for the probability current as the one usually given in quantum mechanics textbooks. We have

$\Psi^{*}\frac{\hat{\vec{p}}}{m} \Psi$

$= - \frac{i \hbar}{m} \Psi^{*} \nabla \Psi$

$= - \frac{i \hbar}{2m} \Psi^{*} \nabla \Psi - \frac{i \hbar}{2m} \Psi^{*} \nabla \Psi$

We now note that since the probability current density must be a real vector, the last two terms above must be real. Therefore they are not affected in any way if we take their complex conjugate. Taking the complex conjugate of the second term in the last equality we get

$\Psi^{*}\frac{\hat{\vec{p}}}{m} \Psi$

$= - \frac{i \hbar}{2m} \Psi^{*} \nabla \Psi + \frac{i \hbar}{2m} \Psi \nabla \Psi^{*}$

$= \frac{i \hbar}{2m}(\Psi \nabla \Psi^{*} - \Psi^{*} \nabla \Psi)$

$= \vec{j}$

This is exactly the expression for the probability current density that appears in textbooks, but rather than introducing it out of nowhere’ at the beginning, we have obtained it naturally as a result of a fluid-mechanical model.

Derivation of Euler’s equations of motion for a perfect fluid from Newton’s second law

Having read a number of highly technical derivations of Euler’s equations of motion for a perfect fluid I feel that the mathematical meanderings tend to obscure the underlying physics. In this note I want to explore the derivation from a more physically intuitive point of view. The dynamics of a fluid element of mass $m$ are governed by Newton’s second law which says that the vector sum of the forces acting on the fluid element is equal to the mass of the element times its acceleration. Thus,

$\vec{F} = m \vec{a}$

The net force $\vec{F}$ can be decomposed into two distinct types of forces, so-called body forces $\vec{W}$ that act on the entire fluid element (e.g., the fluid element’s weight due to gravity) and stresses $\vec{S}$ such as pressures and shears that act upon the surfaces enclosing the fluid element. For the purposes of deriving the differential form of Euler’s equation we will focus on the net force per unit volume acting on the fluid element, $\vec{f}$, which we will decompose into a weight per unit volume $\vec{w}$ and a net stress force per unit volume $\vec{s}$. The weight per unit volume is simply obtained as

$\vec{w} = \rho \vec{g}$

where

$\rho = \frac{m}{V}$

is the mass density of the fluid (i.e., mass per unit volume) and $\vec{g}$ is the acceleration due to gravity. In index notation, the equation for the $i$-th component of the weight per unit volume is

$w^i = \rho g^i$

The net stress force per unit volume, $\vec{s}$, is a little more complicated to derive since it involves the rank-2 stress tensor $\tau^{ij}$. This tensor contains nine components and is usually represented as a $3 \times 3$ symmetric matrix. In Cartesian coordinates the components along the main diagonal, namely $\tau^{xx}$, $\tau^{yy}$ and $\tau^{zz}$, represent normal stresses, i.e., forces per unit area acting orthogonally to the planes whose normal vector is identified by the first superscript, as indicated in the diagram below. (Note that a stress is a force per unit area, so to convert a stress tensor component $\tau^{ij}$ into a force it would be necessary to multiply it by the area over which it acts).

In the diagram each normal stress is shown as a tension, i.e., a normal stress pointing away from the surface. When a normal stress points towards the surface it acts upon, it is called a pressure.

The off-diagonal components of the stress tensor represent shear stresses, i.e., forces per unit area that point along the sides of the fluid element, parallel to these sides rather than normal to them. These shear stresses are shown in the following diagram.

Shear stresses only arise when there is some kind of friction in the fluid. A perfect fluid is friction-free so there are no shear stresses. Euler’s equation only applies to perfect fluids so for the derivation of the equation we can ignore the off-diagonal components of the stress tensor.

The normal stresses along the main diagonal are usually written as

$\tau^{xx} = - p^x$

$\tau^{yy} = - p^y$

$\tau^{zz} = - p^z$

where $p$ stands for pressure and the negative sign reflects the fact that a pressure points in the opposite direction to a tension.

In a perfect fluid the pressure is isotropic, i.e., the same in all directions, so we have

$p^x = p^y = p^z = p$

Therefore the stress tensor of a perfect fluid with isotropic pressure reduces to

$\tau^{ij} = -p \delta^{ij}$

where $\delta^{ij}$ is the Kronecker delta (and may be thought of here as the metric tensor of Cartesian 3-space).

Now suppose we consider the net stress (force per unit area) in the y-direction of an infinitesimal volume element.

The stress on the right-hand face can be approximated using a Taylor series expansion as being equal to the stress on the left plus a differential adjustment based on its gradient and the length $dy$. If we take the stress on the right to be pointing in the positive direction and the one on the left as pointing in the negative (opposite) direction, the net stress in the y-direction is given by

$\tau^{yy} + \frac{\partial \tau^{yy}}{\partial y} dy - \tau^{yy} = \frac{\partial \tau^{yy}}{\partial y} dy$

Similarly, the net stresses in the $x$ and $z$-directions are

$\frac{\partial \tau^{xx}}{\partial x} dx$

and

$\frac{\partial \tau^{zz}}{\partial z} dz$

To convert these net stresses to net forces we multiply each one by the area on which it acts. Thus, the net forces on the fluid element (in vector form) are

$\big(\frac{\partial \tau^{xx}}{\partial x} dxdydz\big) \vec{i}$

$\big(\frac{\partial \tau^{yy}}{\partial y} dxdydz\big) \vec{j}$

$\big(\frac{\partial \tau^{zz}}{\partial z} dxdydz\big) \vec{k}$

The total net force on the fluid element is then

$\big(\frac{\partial \tau^{xx}}{\partial x} \ \vec{i} + \frac{\partial \tau^{yy}}{\partial y} \ \vec{j} + \frac{\partial \tau^{zz}}{\partial z} \ \vec{k}\big) dxdydz$

Switching from tensions to pressures using $\tau^{ij} = -p \delta^{ij}$ and dividing through by the volume $dxdydz$ we finally get the net stress force per unit volume to be

$\vec{s} = -\big(\frac{\partial p}{\partial x} \ \vec{i} + \frac{\partial p}{\partial y} \ \vec{j} + \frac{\partial p}{\partial z} \ \vec{k}\big)$

In index notation, the equation for the $i$-th component of this net pressure per unit volume is written as

$s^i = -\frac{\partial p}{\partial x^i}$

We have now completed the analysis of the net force on the left-hand side of Newton’s second law.

On the right-hand side of Newton’s second law we have mass times acceleration, where acceleration is the change in velocity with time. To obtain an expression for this we observe that the velocity of a fluid element may change for two different reasons. First, the velocity field may vary over time at each point in space. Second, the velocity may vary from point to point in space (at any given time). Thus, we consider the velocity field to be a function of the time coordinate as well as the three spatial coordinates, so

$\vec{v} = \vec{v}(t, x, y, z) = v^x(t, x, y, z) \ \vec{i} + v^y(t, x, y, z) \ \vec{j} + v^z(t, x, y, z) \ \vec{k}$

Considering the $i$-th component of this velocity field, the total differential is

$dv^i = \frac{\partial v^i}{\partial t} \ dt + \frac{\partial v^i}{\partial x} \ dx + \frac{\partial v^i}{\partial y} \ dy + \frac{\partial v^i}{\partial z} \ dz$

so the total derivative with respect to time is

$\frac{dv^i}{dt} = \frac{\partial v^i}{\partial t} \ dt + v^x \frac{\partial v^i}{\partial x} + v^y \frac{\partial v^i}{\partial y} + v^z \frac{\partial v^i}{\partial z}$

where I have used

$v^x = \frac{dx}{dt}$

$v^y = \frac{dy}{dt}$

$v^z = \frac{dz}{dt}$

We can write this more compactly using the Einstein summation convention as

$\frac{dv^i}{dt} = \frac{\partial v^i}{\partial t} + v^j \frac{\partial v^i}{\partial x^j}$

This is then the $i$-th component of the acceleration vector on the right-hand side of Newton’s second law. In component form, therefore, we can write mass times acceleration per unit volume for the fluid element as

$\rho \big(\frac{\partial v^i}{\partial t} + v^j \frac{\partial v^i}{\partial x^j}\big)$

This completes the analysis of the mass times acceleration term on the right-hand side of Newton’s second law.

In per-unit-volume form, Newton’s second law for a fluid element is

$\vec{w} + \vec{s} = \rho \vec{a}$

and writing this in the component forms derived above we get the standard form of Euler’s equations of motion for a perfect fluid:

$\rho g^i - \frac{\partial p}{\partial x^i} = \rho \big(\frac{\partial v^i}{\partial t} + v^j \frac{\partial v^i}{\partial x^j}\big)$