5. The Navier-Stokes equations

Our adventures in fluid dynamics thus far have mostly consisted thus far of considering the ‘macroscopic’ effects of fluids and their interactions with the outside world. We have spent considerable effort over the last few sections in examining the control volume approach, often combined with simplifying assumptions to make the problems more tractable.

Although these are powerful concepts which allow us to leverage this integral approach for even quite complex engineering applications, ultimately we are missing the fine-scaled details: what if we want to know what the fluid velocity is at every point in the domain of interest, for example? To achieve such detail, we need to consider an alternative differential approach, which examines vanishingly small control volumes, and leverages calculus to derive equations that govern fluid flows. If we are able to solve these equations then we could theoretically.

As we shall see, in practice these equation will turn out to be rather complex, and solving them by hand for realistic problems is essentially impossible. Although this may therefore seem like a fruitless task, it gives rise to other, equally powerful strategies. With the equations in hand, we can turn to computational modelling to discretise the equations and solve them via computer. This is the basis for computational fluid dynamics (CFD), which is the subject of the first coursework, and a tool used extensively in modern engineering practice.

5.1. Lagrangian and Eulerian reference frames

To apply a differential approach, we need to consider an infintessimally small ‘blob’ of fluid, which we will say has volume \( dx\, dy\, dz\). The fluid volume will move along with the fluid’s velocity and, as we disussed in previous sections, it will deform along according to that velocity.

All of our concepts of conservation and physical laws apply to this blob of fluid, but we need to do this in the frame of reference of the blob. We call this the Lagrangian frame of reference, which is conceptually clear, but in practice very difficult to use for analysis. The corresponding frame of reference of a passive observer is the Eulerian frame of reference.

A useful analogy for the difference between these frames of reference is a river: the Eulerian approach is akin to sitting by the bank of the river and observing fluid passing through a given point in space, whilst the Lagrangian approach instead is to be in a boat that drifts according to the river’s current.

Naturally, we are typically more interested in the Eulerian frame of reference. However, as with our previous sections, our issue is that all of our notions of conservation and understanding of physical laws of motion revolve around the Lagrangian frame. So we therefore need to figure out how to translate between these frames of reference.

5.1.1. The material derivative

Consider an arbitrary quantity \(B\), which we could use to denote temperature, a component of velocity, or really any desirable quantity. Then we could view \(B\) in either:

• the Eulerian sense, so that \(B\) is a function of space and time in a given domain: i.e. we write \(B(\vec{x}, t)\);

• the Lagrangian sense, so that \(B\) is a function of space and time, but along a given path \(\vec{x}(t) = (x(t), y(t), z(t))\), so that we write \(B(\vec{x}(t), t) = B(x(t), y(t), z(t), t)\).

Notice that in the Lagrangian setting, \(B\) is essentially a function of only \(t\). We can therefore compute its ordinary derivative, and use the chain rule to evaluate this in terms of the derivative of \(\vec{x}\) as

\[\begin{split} \begin{align*} \frac{dB}{dt} &= \frac{\partial B}{\partial t} \frac{dt}{dt} + \frac{\partial B}{\partial x} \frac{dx}{dt} + \frac{\partial B}{\partial y} \frac{dy}{dt} + \frac{\partial B}{\partial z} \frac{dz}{dt} \\ &= \frac{\partial B}{\partial t} + \frac{\partial B}{\partial x} \frac{dx}{dt} + \frac{\partial B}{\partial y} \frac{dy}{dt} + \frac{\partial B}{\partial z} \frac{dz}{dt} \\ &= \frac{\partial B}{\partial t} + \frac{d\vec{x}}{dt}\cdot \nabla B \end{align*} \end{split}\]

where \(\nabla B = (\partial_x B, \partial_y B, \partial_z B)\) is the gradient of \(B\). Notice that, very subtly, the terms on the right now express \(B\) in the Eulerian sense.

Moreover, if we select a path that corresponds to the fluid velocity because we want to follow the trajectory of fluid volumes according to that velocity, then we see that

\[ \frac{d\vec{x}}{dt} = \vec{u} \Rightarrow \boxed{\frac{dB}{dt} = \frac{\partial B}{\partial t} + \vec{u}\cdot\nabla B} \]

The term above is called the material derivative. It is usually denoted using the adjusted term \(D/Dt\). Notice that we can ‘factorise’ the operations above, which gives us the definition:

\[ \boxed{\frac{D}{Dt} = \frac{\partial}{\partial t} + \vec{u}\cdot\nabla} \]

A video example

If you’re looking for clarification on the above, including on Eulerian and Lagrangian frames of reference, I highly recommend the following short video which does a great job of explaining this!

5.1.2. Fluid acceleration

Now we have some connection between the Eulerian and Lagrangian frames through the use of the material derivative, we can use this in the setting of our physical laws. For applying conservation of momentum, we are going to need to describe the acceleration of the fluid.

Using our definition of the material derivative, this is easy: we can just substitute in a component of velocity. For example the \(x\)-component of acceleration is given by:

\[ a_x = \frac{Du}{Dt} = \frac{\partial u}{\partial t} + \vec{u}\cdot \nabla u = \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z} \]

We typically write this in a shorthand manner by collecting together all of the components into a single vector equation:

\[ \vec{a} = \frac{D\vec{u}}{dt} = \frac{\partial \vec{u}}{\partial t} + \vec{u}\cdot\vec{\nabla}\vec{u}, \]

where the vector gradient operator \(\vec{\nabla}\) is given by

\[ \vec{\nabla} = \frac{\partial}{\partial x}\hat{i} + \frac{\partial}{\partial y}\hat{j} + \frac{\partial}{\partial z}\hat{k} \]

Typically we will simply write \(\nabla\) in place of \(\vec{\nabla}\).

Example: computing acceleration

Given an Eulerian velocity field \(\vec{u}(x,y,z,t) = 3t\hat{i} + xz\vec{j} + ty^2\vec{k}\), let’s calculate the total acceleration of a fluid particle.

We essentially break this down into two terms. Firstly we compute the local acceleration term \(\partial_t\vec{u}\) as

\[ \frac{\partial\vec{u}}{\partial t} = \hat{i} \frac{\partial}{\partial x} (3t) + \hat{j} \frac{\partial}{\partial y} (xz) + \hat{k} \frac{\partial}{\partial y} (ty^2) = 3\hat{i} + y^2\hat{k}. \]

Then we can proceed to compute the second part of the acceleration, which we call the convective acceleration:

\[\begin{split} \begin{align*} u \frac{\partial\vec{u}}{\partial x} &= 3t\frac{\partial}{\partial x} (3t\hat{i} + xz\vec{j} + ty^2\vec{k}) = 3txz\vec{j} \\ v \frac{\partial\vec{u}}{\partial y} &= xz\frac{\partial}{\partial y} (3t\hat{i} + xz\vec{j} + ty^2\vec{k}) = 2txyz\vec{k} \\ w \frac{\partial\vec{u}}{\partial z} &= ty^2\frac{\partial}{\partial z} (3t\hat{i} + xz\vec{j} + ty^2\vec{k}) = txy^2\vec{j} \\ \end{align*} \end{split}\]

Thus,

\[ \vec{a} = \frac{D\vec{u}}{Dt} = (3\hat{i} + y^2\hat{k}) + (3txz\vec{j}) + (2txyz\vec{k}) + (txy^2\vec{j}) = 3\vec{i} + (3tz + txy^2)\vec{j} + (y^2 + 2txyz)\vec{k}. \]

5.2. Conservation of mass

Now we can begin to consider our conservation laws in this differential setting. Consider that we have an infinitessimal control volume of dimensions \(dx \times dy \times dz\) which is fixed in space, as shown in the figure below.

Fig. 5.1 A differential volume.

Since the volume is very small, we can consider the flow through each face to be uniform. Thus, if we now take our volume to be a control volume, we can apply some of our integral relations from the previous sections. In particular then we have that, for a general fluid which has potential variable density, conservation of mass on the volume is given by the continuity equation

\[ \int_{CV} \frac{\partial\rho}{\partial t} dV + \sum_i (\rho_i A_i V_i)_{\mathrm{out}} - \sum_i(\rho_i A_i V_i)_{\mathrm{in}} = 0. \]

Now, if the volume is small then our continuum hypothesis can apply, and we say that the material properties on the control volume are constant. This means that

\[ \int_{CV} \frac{\partial\rho}{\partial t} dV = \frac{\partial\rho}{\partial t} dV = \frac{\partial\rho}{\partial t} dx\,dy\,dz. \]

Moreover, our cube has 3 inlets and 3 outlets. In the figure above, we outline the mass flow through the faces in the \(x\) direction. Both faces have an area of \(dy\, dz\); the inlet mass flow is thus \(\rho u dy\, dz\) and the outlet mass flow is \(\left[\rho u + dx\,\partial_x (\rho u)\right] dy\, dz\).

If we do this for each of the faces, then substitute this into our conservation of mass equation above, we have

\[ \frac{\partial\rho}{\partial t} dx\,dy\,dz + \frac{\partial}{\partial x}(\rho u) dx\,dy\,dz + \frac{\partial}{\partial y}(\rho v) dx\,dy\,dz + \frac{\partial}{\partial z}(\rho w) dx\,dy\,dz = 0 \]

which simplifies immediately to

\[ \boxed{\frac{\partial\rho}{\partial t}+ \frac{\partial}{\partial x}(\rho u) + \frac{\partial}{\partial y}(\rho v) + \frac{\partial}{\partial z}(\rho w) = 0.} \]

This is the continuity equation in differential form. More compactly, we could write this using vector notation as

\[ \frac{\partial\rho}{\partial t} + \nabla\cdot(\rho\vec{u}) = 0. \]

where \(\nabla \cdot \vec{u} = \partial_x u + \partial_y v + \partial_z w\) is the divergence of \(\vec{u}\).

5.2.1. Incompressible flows

Where the flow is incompressible, we know that the density \(\rho\) is constant. Then \(\partial_t \rho = 0\) and we can rewrite the continuity equation as

\[ \nabla\cdot(\rho u) = 0 \Rightarrow \rho \nabla\cdot\vec{u} = 0 \Rightarrow \boxed{\nabla\cdot\vec{u} = 0}. \]

We say then in this case that the velocity field \(\vec{u}\) is divergence-free.

Calculating the form of a velocity field

An incompressible velocity field \(\vec{u} = (u,v,w)\) is known to have components \(u = a(x^2-y^2)\) and \(w=b\), where \(a\) and \(b\) are constants. What form must the component \(v\) be?

We know that the velocity field is divergence free, so

\[ \nabla\cdot\vec{u} = 2ax + \frac{\partial v}{\partial y} = 0 \Rightarrow \frac{\partial v}{\partial y} = -2ax. \]

Therefore, we can integrate with respect to \(y\) to yield

\[ v = -2axy + f(x,z,t) \]

where \(f(x,z,t)\) is an unknown function that depends on all variables except \(y\).

5.3. Conservation of (linear) momentum

Let’s now turn to the conservation of momentum, and use the same approach as in the previous section as we did for mass. We have from our previous discussion of conservation of momentum on a control volume that the net force \(\vec{F}\) must be balanced by the change in momentum flux: i.e. Newton’s second law. Using the same approach we did in the previous section, we can therefore write that

\[ \vec{F} = \frac{\partial}{\partial t}\left(\int_{\mathrm{CV}} \vec{u}\rho \,dV\right) + \sum_i (\dot{m}_i \vec{u}_i)_{\mathrm{out}} - \sum_i(\dot{m}_i \vec{u}_i)_{\mathrm{in}} \]

where \(\dot{m}\) is the mass flow through an inlet or an outlet.

5.3.1. Momentum fluxes

Let’s first ignore the forces, and consider only things on the right hand side. We make the same assumptions about material properties being essentially constant on the control volume due to its small size, and we can compute the inlet and outlet momentum fluxes according to the table below:

Face	Inlet momentum flux	Outlet momentum flux
\(x\)	\(\rho u\vec{u}dy\,dz\)	\(\left[\rho u\vec{u} + \frac{\partial}{\partial x}(\rho u\vec{u})\right]dy\,dz\)
\(y\)	\(\rho v\vec{u}dx\,dz\)	\(\left[\rho v\vec{u} + \frac{\partial}{\partial y}(\rho v\vec{u})\right]dx\,dz\)
\(z\)	\(\rho w\vec{u}dx\,dy\)	\(\left[\rho w\vec{u} + \frac{\partial}{\partial z}(\rho w\vec{u})\right]dx\,dy\)

Substituting all of this into the equation above, and again cancelling \(dx\, dy\, dz\) we obtain the expression

\[ \frac{\partial}{\partial t} (\rho\vec{u}) + \frac{\partial}{\partial x} (\rho u\vec{u}) + \frac{\partial}{\partial y} (\rho v\vec{u}) + \frac{\partial}{\partial z} (\rho w\vec{u}). \]

Now, we can essentially apply the product rule for each of the partial derivatives in terms of \(\rho\) and \(\vec{u}\) to obtain the expression

\[ \vec{u}\left[ \frac{\partial \vec{u}}{\partial t} + \nabla\cdot(\rho \vec{u}) \right] + \rho \left[\frac{\partial \vec{u}}{\partial t} + \vec{u}\cdot\nabla \vec{u}\right] = \underbrace{\vec{u}\left[ \frac{\partial \vec{u}}{\partial t} + \nabla\cdot(\rho \vec{u}) \right]}_{\mathrm{continuity\ equation}} + \rho\underbrace{\frac{D\vec{u}}{Dt}}_{\mathrm{acceleration}}. \]

Since we know the continuity equation must be zero, this leaves us with the final equation

\[ \vec{F} = \rho dx\,dy\,dz \frac{D\vec{u}}{dt} \]

5.3.2. Forces due to gravity, viscosity and pressure

Phew – unfortunately though, we’re only half way there! Now we need to turn our attention to the net force on the left hand side. We know from our extended discussion so far that this is due to three things: gravity (a body force), viscosity (a shear force) and pressure (a normal force).

5.3.2.1. Graviational force

This is the easiest force to deal with. We know that this is a body force, so just like in the very first section of these notes, we can write that

\[ \vec{F}_g = \rho \vec{g} dx\,dy\,dz \]

Typically we would take the \(z\)-direction as upwards, in which case \(\vec{g} = -g\hat{k}\).

5.3.2.2. Pressure and viscous forces

In our discussion of internal flows, we made lots of assumptions about the flow: like that the flow was laminar and had a velocity profile that was one-dimensional, so that shear forces due to viscosity only operated on, for example, the top and bottom of a small volume.

Unfortunately, in reality, every face of the volume is going to have a stress imposed on it due to both viscosity shear and pressure. In total there will be nine of them: stresses in the \(x\), \(y\), and \(z\) directions that act normal to the \(x\), \(y\) and \(z\) directions. We write each of these as a variable \(\sigma_{ij}\), which denotes the stress in direction \(j\) on a face normal to the \(i\) axis. This is visualised in the following figure.

Fig. 5.2 The stress tensor \(\sigma_{ij}\).

We write these stresses not as individual components (typically), but instead as a matrix so that

\[\begin{split} \sigma_{ij} = \begin{bmatrix} \tau_{xx} - p & \tau_{yx} & \tau_{zx} \\ \tau_{xy} & \tau_{yy} - p & \tau_{zy} \\ \tau_{xz} & \tau_{yz} & \tau_{zz} - p \end{bmatrix} \end{split}\]

This is called a tensor, since we rely not only on a single index to define a component, as in a vector, but instead two components. We shall not really consider tensor analysis any further in this course, but it is useful for you to be exposed to the concepts.

The \(\sigma_{ij}\) stress tensor consists of two components:

• the pressure force which, being a normal force that is equal in all directions simply acts on the \(x\), \(y\) and \(z\) faces in the same directions;

• the \(\tau_{ij}\) terms, which denote the shear forces acting on each face in each direction.

These stresses will eventually give rise to forces, which we can compute since we know the area of the faces in each direction. Notably, applying the same logic as before, we see that the change in stresses are what is responsible for forces on the control volume. Considering only the \(x\)-direction (for clarity), we arrive at the following figure:

Fig. 5.3 The forces acting on the control volume due to the stress tensor \(\sigma_{ij}\) in the \(x\)-direction only.

Calculating the surface force in the \(x\)-direction by combining these terms therefore yields:

\[ \vec{F}_{\mathrm{surf},x} = \left[ -\frac{\partial p}{\partial x} + \frac{\partial}{\partial x}(\tau_{xx}) + \frac{\partial}{\partial x}(\tau_{yx}) + \frac{\partial}{\partial x}(\tau_{zx})\right] dx\,dy\,dz. \]

The remaining \(y\) and \(z\)-directions follow in exactly the same fashion (albeit with a tedious amount of calculation), giving

\[\begin{split} \begin{align*} \vec{F}_{\mathrm{surf},y} &= \left[ -\frac{\partial p}{\partial y} + \frac{\partial}{\partial x}(\tau_{xy}) + \frac{\partial}{\partial x}(\tau_{yy}) + \frac{\partial}{\partial x}(\tau_{zy})\right] dx\,dy\,dz\\ \vec{F}_{\mathrm{surf},z} &= \left[ -\frac{\partial p}{\partial z} + \frac{\partial}{\partial x}(\tau_{xz}) + \frac{\partial}{\partial x}(\tau_{yz}) + \frac{\partial}{\partial x}(\tau_{zz})\right] dx\,dy\,dz\\ \end{align*} \end{split}\]

Mercifully, we are finished. We can return to our original conservation equations and substitute in all of these forces:

\[ \vec{F}_g + \vec{F}_{\mathrm{surf}} = \rho dx\,dy\,dz \frac{D\vec{u}}{dt} \]

which, in the fully expanded form, gives the incredibly complex equations:

\[\begin{split} \begin{align*} \rho g_x - \frac{\partial p}{\partial x} + \frac{\partial \tau_{xx}}{\partial x} + \frac{\partial \tau_{yx}}{\partial y} + \frac{\partial \tau_{zx}}{\partial z} &= \rho\left(u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} + w \frac{\partial u}{\partial z} \right) \\ \rho g_z - \frac{\partial p}{\partial y} + \frac{\partial \tau_{xy}}{\partial x} + \frac{\partial \tau_{yy}}{\partial y} + \frac{\partial \tau_{zy}}{\partial z} &= \rho\left(u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} + w \frac{\partial v}{\partial z} \right) \\ \rho g_z - \frac{\partial p}{\partial z} + \frac{\partial \tau_{xz}}{\partial x} + \frac{\partial \tau_{yz}}{\partial y} + \frac{\partial \tau_{zz}}{\partial z} &= \rho\left(u \frac{\partial w}{\partial x} + v \frac{\partial w}{\partial y} + w \frac{\partial w}{\partial z} \right) \end{align*} \end{split}\]

Typically, we would not write the equations in such a manner. If we define the divergence of \(\boldsymbol\tau = \tau_{ij}\) as

\[ \nabla\cdot\boldsymbol{\tau} = \hat{i} \left(\frac{\partial\tau_{xx}}{\partial x} + \frac{\partial\tau_{yx}}{\partial y} + \frac{\partial\tau_{zx}}{\partial z} \right) + \hat{j} \left(\frac{\partial\tau_{xy}}{\partial x} + \frac{\partial\tau_{yy}}{\partial y} + \frac{\partial\tau_{zy}}{\partial z} \right) + \hat{k} \left(\frac{\partial\tau_{xz}}{\partial x} + \frac{\partial\tau_{yz}}{\partial y} + \frac{\partial\tau_{zz}}{\partial z} \right) \]

Then this simplifies to

\[ \boxed{\rho\vec{g} - \nabla p + \nabla\cdot\boldsymbol{\tau} = \rho\frac{D\vec{u}}{dt}} \]

which looks far nicer, albeit hiding a lot of complexity under the notation! These then are the full equations of state for any fluid, compressible or incompressible.

Warning

Don’t forget that we actually have to make sure that this is combined also with the continuity equation!

5.3.3. Incompressible Newtonian flows

Since the equation above is, frankly, completely horrendous to deal with in any reasonable manner, let’s bring the level of complexity down to ‘moderately horrendous’ by considering a simplification: incompressible Newtonian fluids.

In this case, we know that the shear stresses owing to viscosity are proportional to the rate of shear strain. In one dimension we could write this approximately as

\[ \tau = \mu\frac{du}{dx} \]

where \(\mu\) is the dynamic viscosity. In a general three-dimensional setting, we actually see pretty much the same relationship. The slightly complex way to write this is that \(\boldsymbol\tau = \mu(\nabla\vec{u} + \nabla\vec{u}^{T})\). This gives three diagonal components of \(\boldsymbol\tau\),

\[ \tau_{xx} = 2\mu\frac{du}{dx}, \quad \tau_{yy} = 2\mu\frac{dv}{dy}, \quad \tau_{zz} = 2\mu\frac{dw}{dz} \]

and three unique off-diagonal components owing to the symmetry of the expression:

\[ \tau_{xy} = \tau_{yx} = \mu\left(\frac{du}{dy} + \frac{dv}{dx} \right), \quad \tau_{xz} = \tau_{zx} = \mu\left(\frac{dw}{dx} + \frac{du}{dz} \right), \quad \tau_{xz} = \tau_{zx} = \mu\left(\frac{dv}{dz} + \frac{dw}{dy}. \right) \]

If we plug these into the equations above, and divide everything by the now-constant density \(\rho\), we are left with the Navier-Stokes equations

\[ \boxed{\frac{\partial\vec{u}}{\partial t} + \vec{u}\cdot\nabla\vec{u} = -\frac{1}{\rho} \nabla p + \nu\nabla^2\vec{u} + \vec{g}} \]

These are still extremely difficult to solve, in general, but they are at least a little more tractable to deal with! As we mentioned above, computational fluid dynamics software discretises these equations of state, and transforms them into essentially linear algebra problems. This is the subject of your first coursework!

Warning

You should note that you still need to solve the continuity equation. Looking at the unknowns, we see that we have four: \(u\), \(v\), \(w\) and \(p\). However, the conservation of momentum above only gives three equations: without the conervation of mass, we won’t be able to solve these problems in general.

5.4. Laminar flow between two plates

After all of this effort, we should probably confirm that our equations validate what we already know to be true! For one last time, let us consider laminar flow between two parallel plates.

Example: Navier-Stokes for flow between two flat plates

To apply Navier-Stokes, we need to again leverage the assumption that our velocity vector is only unidirectional and is only a function of \(y\), as shown in the figure below, so that \(\vec{u} = u(y)\hat{i}\). We also assume it is steady.

Fig. 5.4 A channel flow again, with walls at \(y=-L, L\).

In this case, the divergence of \(u\) is automatically zero: so the continuity equation doesn’t really tell us anything. Let’s instead then consider the conservation of momentum. Since \(v\) and \(w\) are both zero, we only consider the first equation in terms of \(y\). Writing this out, we have

\[ \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z} = -\frac{1}{\rho}\frac{\partial p}{\partial x} + \nabla\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} \right) \]

Now, most of these terms are zero! The the first term is zero since the flow is steady; all of the convective terms are zero since \(v=w=0\) and \(\partial_x u = 0\); and all but one of the Laplacian terms is zero. This then gives that

\[ \frac{1}{\rho} \frac{\partial p}{\partial x} = \nu\frac{\partial^2 u}{\partial y^2} \Rightarrow \frac{\partial p}{\partial x} = \mu\frac{\partial^2 u}{\partial y^2} \]

Since \(p\) will only be a function of \(x\), and \(u\) only a function of \(y\), these could be expressed instead as ordinary derivatives, i.e.

\[ \frac{dp}{dx} = \mu\frac{d^2 u}{dy^2} \]

Looking back at the section on internal flows, this is precisely the expression we derived via control volume analysis. So we won’t repeat this analysis as we’ll only get the same conclusions!

5.5. Summary & further reading

This section has shown you:

• the Lagragian and Eulerian views of a fluid;

• how to calculate material derivatives;

• the continuity of mass in a differential setting;

• the continuity of momentum in a differential setting;

• and the Navier-Stokes equations.

This is the final section on fluids! Our next section will introduce thermodynamics.

Further reading

For further reading, consult chapter 4 of White. In particular:

• section 4.1 for discussion of the acceleration of a fluid and the Eulerian/Lagrangian approaches;

• section 4.2 for differential forms of the continuity of mass;

• and section 4.3 for the differential form the continuity of momentum and Navier-Stokes equations.