Surprisingly, most people don’t, and by people I mean those with either a degree or relevant job experience in the area of computational fluid dynamics (CFD). Did you also know that there was a third French scientist (whose name is almost lost to history) who actually contributed much more to the development of the Navier-Stokes equations than Navier did? In this article, I am going through a step by step derivation of the Navier-Stokes equation at a level of detail you will not have seen before, and I have included historic annotation to understand how the equation came about. More importantly, I explain exactly what hypothesis Navier and Stokes had to make the momentum equation account for viscous forces and how their molecular viewpoints translate into partial differential equations. By researching for this article I could not find a single source explaining the actual hypothesis of Navier and Stokes. It is usually just given in its final form without explanation. You will find here a detailed step-by-step derivation of the Navier-Stokes equations and you will fully understand each term by the end of this page. Let’s get started then.
Table of Contents
Deriving Reynolds Transport Theorem
Before we can derive the full Navier-Stokes equations, we need to first investigate how a flow property evolves over time. For this, we consider a volume V_0 which is depicted for some time t below
At time t, we have a control volume of some shape which evolves over time. Let’s assume that for some time y we define V(t)=V_0 and, after some time t+\Delta t, the control volume changes to V_0+\Delta V. Thus, we simply have provided a way to describe mathematically that a certain (control) volume can change over time. You can think of this control volume as a collection of atoms (or the convex hull of a group of atoms). After some time t+\Delta t, the atoms have moved, each atom with its own velocity vector and direction. Thus, the control volume (or convex hull) must have changed as a result. This is exactly what we see depicted above.
Consider now an extensive flow property (a property that depends on the system size such as mass and volume) \Phi, which is conserved within the control volume. We can state that it changes according to
\Phi = \int_V \phi \mathrm{d}V.This makes intuitive sense. Extensive properties depend on the system size, e.g. the mass of a 1 m^3 domain will be different than that of a 10 m^3 domain, assuming the same density. Thus, the above equation simply states that when we increase or decrease the volume of the system through the integral \int_V\mathrm{d}V, then this change in volume is directly proportional to changes in \Phi. There is a proportionality factor \phi, which depending on \Phi will take different forms.
If we consider \Phi to be mass, then \phi becomes the density of the fluid referred to as \rho. This is shown below
m = \int_V \rho\cdot \mathrm{d}V.This makes sense as a density multiplied by a volume provides a mass. This can be also easily seen by just looking at the units
\left[kg\right]=\left[\frac{kg}{\sout{m^3}}\right]\cdot\left[\sout{m^3}\right]Here, we are assuming that \Phi is our mass in kilograms and \phi our density in kilograms per cubic meter. The beauty of the above equation for \Phi is that it is general so we can apply it to other quantities as well.
If, for example, we now ask \Phi to represent the impulse m\mathbf{u}, then \phi becomes the impulse density \rho\mathbf{u}. This is shown in the equation below
m\mathbf{u}=\int_V \rho\mathbf{u}\cdot\mathrm{d}V.We can show, again, by just looking at units that this is correct
\left[kg\cdot\frac{m}{s}\right]=\left[\frac{kg}{\sout{m^3}}\cdot\frac{m}{s}\right]\cdot\left[\sout{m^3}\right]Alternatively, we can use the definition we saw before where \Phi represented the mass m and simply multiply this definition by the velocity vector \mathbf{u}. Both of these definitions will become important once we start deriving our conservation laws of fluid dynamics.
Next, we are interested in how \Phi changes over time. This will allow us to monitor \Phi over time and as the control volume V_0 is changing. To do so, we simply add a time derivative on either side to arrive at
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}t}\int_{V(t)}\phi\mathrm{d}V.This equation does provide us now with a time history, however, we have no idea (or means to measure) how the control volume is changing over time and thus we can’t solve the integral on the right-hand side, as it depends on V(t). Whenever we are faced with a situation where we don’t know how a quantity is changing over time or space, we can create a Taylor series for the quantity of interest around some finite point in time or space and some small increment around that point in space or time. For example, a simple Taylor series for \phi in space would give us
\phi(x+\Delta x) = \phi(x) + \frac{\partial \phi}{\partial x}\frac{\Delta x}{1!} + \mathcal{O}(\Delta x)^2.In time, we could develop a Taylor series around \phi in the following way
\phi(t+\Delta t) = \phi(t) + \frac{\partial \phi}{\partial t}\frac{\Delta t}{1!} + \mathcal{O}(\Delta t)^2.Both Taylor series will provide us with knowledge of how \phi will change close to x or t. We can rearrange either one of these Taylor series and write (here, using the Taylor series in time and noting that 1!=1)
\frac{\phi(t+\Delta t) - \phi(t)}{\Delta t} = \frac{\partial \phi}{\partial t} + \mathcal{O}(\Delta t)^2.Ignoring the higher order terms, which we have lumped together in \mathcal{O}(\Delta t)^2, we can further write the following approximation
\frac{\partial \phi}{\partial t} \approx \frac{\phi(t+\Delta t) - \phi(t)}{\Delta t}This approximation is valid for small \Delta t values. If we wanted to have an exact derivative, then we require that \Delta t\rightarrow 0, which we can express in the form of a limit, i.e. \lim\limits_{\Delta t\rightarrow 0}.
Returning now to our time-dependent formulation of \Phi, i.e.
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}t}\int_{V(t)}\phi\mathrm{d}V,we can now make use of the Taylor series expansion in time we derived above and instead of looking at the volume change over some arbitrary interval \Delta t, we check how the volume changes in time, i.e. V(t) as we approach the limit of zero. This can be expressed as
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\lim\limits_{\Delta t\rightarrow 0}\frac{1}{\Delta t}\left\{ \int_{V_0 + \Delta V} \phi(t + \Delta t)\mathrm{d}V - \int_{V_0} \phi(t)\mathrm{d}V \right\},Here, we have simply introduced the approximation of \partial \phi /\partial t and placed \Delta t outside the integration as it remains constant. Note that because the volume can change for some t + \Delta t, the integration is carried out over V_0 + \Delta V, while the integration at t can be taken over V_0.
We can expand the first of the two integrals and arrive at
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\lim\limits_{\Delta t\rightarrow 0}\frac{1}{\Delta t}\left\{ \int_{V_0} \phi(t + \Delta t)\mathrm{d}V + \int_{\Delta V} \phi(t + \Delta t)\mathrm{d}V - \int_{V_0} \phi(t)\mathrm{d}V \right\}.This allows us to group together integrals of the same bounds (i.e. here \int_{V_0}). By doing so, we can write the above equation as
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\lim\limits_{\Delta t\rightarrow 0}\frac{1}{\Delta t}\left\{ \int_{V_0} \left[ \phi(t + \Delta t) - \phi(t) \right] \mathrm{d}V + \int_{\Delta V} \phi(t + \Delta t)\mathrm{d}V \right\}.Let’s think now about what we want to achieve. We want to derive a conservation law that we can use to derive arbitrary conservation equations, for example for mass and momentum conservation. These types of equations are typically given in the integral or derivative form. Thus, in order to arrive at such an equation, we need to eliminate the limit from our equation. This is what we will focus on next. To get started, we expand the above equation and write the limit in front of each integral, This will allow us to evaluate each integral on its own. We have
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\lim\limits_{\Delta t\rightarrow 0}\frac{1}{\Delta t}\int_{V_0} \left[ \phi(t + \Delta t) - \phi(t) \right] \mathrm{d}V + \lim\limits_{\Delta t\rightarrow 0}\frac{1}{\Delta t} \int_{\Delta V} \phi(t + \Delta t)\mathrm{d}V.Let us now look at the first integral. First, let’s reintroduce \Delta t back into the integral. Remember a constant can be written in front or inside the integral. This results in
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\lim\limits_{\Delta t\rightarrow 0}\int_{V_0} \frac{\phi(t + \Delta t) - \phi(t)}{\Delta t} \mathrm{d}V + \lim\limits_{\Delta t\rightarrow 0}\frac{1}{\Delta t} \int_{\Delta V} \phi(t + \Delta t)\mathrm{d}V.A few equations ago, we have already established the following relationship through the Taylor series
\frac{\partial \phi}{\partial t} = \lim\limits_{\Delta t\rightarrow 0}\frac{\phi(t+\Delta t) - \phi(t)}{\Delta t}Thus, we can simply rewrite the first integral as the time derivative which will remove the limit from the integral. This is given as
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\int_{V_0} \frac{\partial \phi}{\partial t} \mathrm{d}V + \lim\limits_{\Delta t\rightarrow 0}\frac{1}{\Delta t} \int_{\Delta V} \phi(t + \Delta t)\mathrm{d}V.We have successfully removed the limit for the first integral and for the second, we see that what we want to do here is to replace \phi(t + \Delta t) with a Taylor series approximation as we have done before. This should drive home the point that whenever we don’t know how a quantity is changing, use a Taylor series around that point in space or time (or both, Taylor series can be established for multi-variable problems) and then work with the terms you get. In our case, if we insert the Taylor series expansion for \phi(t + \Delta t), we get
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\int_{V_0} \frac{\partial \phi}{\partial t} \mathrm{d}V + \lim\limits_{\Delta t\rightarrow 0}\frac{1}{\Delta t} \int_{\Delta V}\left[ \phi(t) + \frac{\partial \phi}{\partial t}\frac{\Delta t}{1!} + \mathcal{O}(\Delta t)^2 \right]\mathrm{d}V. Now let’s think about the terms involved and what happens as we let \Delta t \rightarrow 0. Any terms involving \Delta t have to become zero by definition, thus the second term, i.e. \frac{\partial \phi}{\partial t}\frac{\Delta t}{1!} becomes zero as \Delta t. Any other higher order term containing \Delta t also has to become zero and thus we can say that \mathcal{O}(\Delta t)^2 also becomes zero. The first term, however, i.e. \phi(t) is not zero and can take any shape or form in time and thus remains in the integral. Thus, we can write the second integral simply as
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\int_{V_0}\frac{\partial \phi}{\partial t}\mathrm{d}V+\lim\limits_{\Delta t\rightarrow 0}\frac{1}{\Delta t}\int_{\Delta V} \phi(t)\mathrm{d}V.We were not able to get rid of the limit entirely; we still have a term of 1/\Delta t. If we allowed the limit of this term to go to zero, then the second integral would vanish. While this is mathematically correct, it is the equivalent of saying 0=0, i.e. it is a trivial solution and does not help us understand how some quantity \Phi changes in time and space (carry out the integration above for the first integral and drop the second and you see that we indeed have found a trivial solution). Therefore, we need to find a way to eliminate the \Delta t first.
One way of doing so is to transform the volume integral \mathrm{d}V to a surface integral \mathrm{d}A. To do so, remember that a volume can be constructed from an area that is swept through space for some distance L, or
V=A\cdot L.
We can always express a distance L as a product of a speed and time interval over which we travel with a certain speed. Knowing that we are interested still in some finite time interval, we can replace L as
V=A\cdot\mathbf{u}\Delta t.The equation still holds true for finite volumes and finite areas, i.e.
\mathrm{d}V=\mathrm{d}A\cdot\mathbf{u}\Delta t.We can insert this definition into our previous integral equation and we obtain
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\int_{V_0}\frac{\partial \phi}{\partial t}\mathrm{d}V+\lim\limits_{\Delta t\rightarrow 0}\frac{\Delta t}{\Delta t}\int_{A}\phi\mathbf{u}\mathrm{d}A,or, realising that \Delta t / \Delta t is unity and thus will no longer influence the limit, we can drop the limit entirely from our equation to arrive at
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\int_{V_0}\frac{\partial \phi}{\partial t}\mathrm{d}V+\int_{A}\phi\mathbf{u}\mathrm{d}A.The above equation is known as the Reynolds Transport Theorem and tells us how a quantity \Phi changes in space and time within a fluid. We can rewrite the Reynolds transport theorem for the second integral using Gauss’ theorem, which allows us to transform surface integrals into volume integrals (and vice versa), which provides us with
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\int_{V_0}\frac{\partial \phi}{\partial t}\mathrm{d}V+\int_{V}\nabla\cdot\phi\mathbf{u}\mathrm{d}V.This form is particularly useful if we want to derive a differential form of the Reynolds transport theorem. This can be achieved by requiring now that the volumes we are looking at are going to zero, i.e. V \rightarrow 0. In this case, we obtain
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\frac{\partial \phi}{\partial t}+\nabla\cdot\phi\mathbf{u}.In the literature, you will find different names for \mathrm{d}\Phi(t)/\mathrm{d}t, it is either referred to as the total derivative, the substantial derivative or, the material derivative. All of these denote the same thing and when you read the literature, you may come across these different names.
This website exists to create a community of like-minded CFD enthusiasts and I’d love to start a discussion with you. If you would like to be part of it, sign up using the link below and you will receive my OpenFOAM quick reference guide, as well as my guide on Tools every CFD developer needs for free.
Join now
Deriving the Conservation Laws of Fluid Dynamics
Now that we have spent some time deriving the Reynolds Transport Theorem, we can put it to use to derive the governing equations of fluid dynamics. In this section, we look at the conservation of mass and momentum, which form the basic building blocks for any CFD simulation.
Conservation of Mass
We have already done all the heavy lifting and derived the Reynolds Transport theorem; deriving a conservation law for mass simply requires applying this theorem.
We have already alluded to before that if want \Phi to be the mass, then we need to set \phi=\rho. Thus, we write the Reynolds Transport Theorem as
\frac{\mathrm{d}m}{\mathrm{d}t}=\frac{\partial \rho}{\partial t}+\nabla\cdot\rho\mathbf{u}=0.The right-hand side is set to zero because that is what a conservation law is telling us. It states that all terms in the conservation law added up should result in a value of zero, in this way all quantities are conserved and no mass, in this case, is produced nor destroyed (this assumption always requires considering some form of a closed system).
Now all conservation laws have one thing in common and that is that the property they describe (here the mass) is conserved. This is stating the obvious – it’s in the name – but this means that the time derivative of the mass has to be zero. Therefore, the conservation of mass, or continuity equation, can be written as
\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf{u})=0.This equation describes, in general, how mass is conserved. We will take special considerations into account at the end of this article to show how it may change for different flow regimes.
For an incompressible fluid, we say that the density is constant and so with \rho=const., we have
\nabla\cdot(\rho\mathbf{u})=0.Since \rho=const., we can also divide by it, which will further simplify our equation to
\nabla\cdot\mathbf{u}=0.This is the mass conservation law for an incompressible fluid and is only applicable to these types of flows. It is an important concept that will help us simplify other conservation laws encountered in incompressible flows.
Conservation of Momentum
The conservation of momentum is derived from Newton’s second law which is given as
\mathbf{F}=m\mathbf{a},where \mathbf{F} is the sum of the internal and external forces acting on the system and \mathbf{a} the acceleration. The question then becomes, which forces are acting on our system? Typically this is separated into two broad categories; internal and external forces and this can be written as
\mathbf{F} = \mathbf{F}_{internal} + \mathbf{F}_{external}.Let’s look at the external forces first then. The external forces describe anything that acts on our system but that is not explicitly captured by our Navier-Stokes equations. It acts on each fluid particle (or, if we run a CFD simulation, on each cell of our computational domain). Think of two immiscible fluids, where the heavy fluid is sitting on top of the light fluid as shown in the video below
But think about where the “force” is coming from in order for the heavier fluid to penetrate into the lighter fluid. If we just solve the Navier-Stokes equations without any external forces, nothing would happen. But at this point, we probably can all agree that it is gravity that is pulling the heavier fluid into the lighter fluid. Thus, gravity is one such external force that we can consider. But this also provides an additional point to consider; if we do not expect gravity to be important for us in our situation (i.e. simulation), then it is OK to neglect the external forces such as gravity. But when we derive the conservation of momentum, we want to keep all those terms and then drop them later if we think this is appropriate.
Similar to gravity, there are other forces acting on the fluid which are typically termed body forces, i.e.
\mathbf{F}_{external} = \mathbf{F}_{gravity} + \mathbf{F}_{body}.Examples of the body forces are magnetic fields which influence the trajectory of magnetised fluid particles or the Coriolis force, which imparts momentum in a specific direction onto the fluid particles. Essentially anything that would translate into additional momentum on fluid particles is lumped together into the body forces (and in some derivation, gravity is part of that body force term).
The internal forces are a bit more complicated and require some further considerations. Let’s think about a simplified system, we have two solid blocks, one sitting on top of the other. The block on top will assert a pressure force on the block below (similar to the heavier and lighter fluid situation above) and the block below will assert the same pressure force back according to Newton’s third law (when two bodies interact with each other, they will apply an equal force to each other with opposite direction). So one such internal force is the pressure (force) itself. Now let’s imagine we were to push that block on top forwards so that we have friction between these two blocks. This friction force is due to the upper block shearing over the bottom block and we feel that as a resistance when pushing the upper block. This is referred to as the shear force or more commonly the shear stresses. This is what forms part of the internal forces and we can say that it depends on forces due to pressure and shear stresses, i.e.
\mathbf{F}_{internal} = \mathbf{F}_{pressure} + \mathbf{F}_{shear stresses}.In summary, we have forces of all different forms and shapes. In the following derivation, we will carry over internal, gravity, and body forces for the remaining derivation. Let’s look at these terms in turn.
The internal force contains the pressure and shear stresses. From a dimensions-based analysis, we know that shear stresses are defined, according to Newton, as
\tau=\mu\frac{\partial\mathbf{u}}{\partial\mathbf{x}}\rightarrow\left[\frac{kg}{m\cdot s}\right]\cdot\left[\frac{m}{m\cdot s}\right]=\left[\frac{kg}{m\cdot s^2}\right]=\left[\frac{kg\cdot m}{m^2\cdot s^2}\right]=\left[\frac{N}{m^2}\right].If we now were to take the derivative of \tau as in \nabla\cdot\tau, then we would obtain units of N/m^3. Integrating that over a volume would result in N, which is the correct unit for a force. Thus, we can write that
\mathbf{F}_{internal} = \int_{V}\nabla\cdot\sigma\mathrm{d}V,where we have chosen to represent all internal forces by \sigma, which includes not just the shear stresses but also the pressure forces. We will return to this point later. We do the volume integral here as the internal force is not a constant but may vary over our entire system’s volume. Integrating over all contributions will result in the overall internal force.
Gravity has units of m/s^2 and thus we need to multiply it with some mass to get units of kg\cdot m/s^2=N. But since we are also going to integrate over the system’s volume (similar to the internal forces), we need to divide this force by units of m^3, which means we need to multiply our gravitational force by something of units kg/m^3, which is just the density \rho. Therefore, the gravitational force becomes
\mathbf{F}_{gravity} = \int_{V}\rho\mathbf{g}\mathrm{d}V.The body forces are provided as forces themselves and have units of N. Thus, we have to simply divide them by the system’s volume to account for the volume integral as in the cases before, and arrive at
\mathbf{F}_{body} = \frac{1}{V}\int_{V}\mathbf{f}_{body}\mathrm{d}V.Thus, in summary, we can write the overall force as a superposition of the internal, gravitational, and body forces as
\mathbf{F} = \mathbf{F}_{internal} + \mathbf{F}_{gravity} + \mathbf{F}_{body} =
\int_{V}\nabla\cdot\sigma\mathrm{d}V + \int_{V}\rho\mathbf{g}\mathrm{d}V +
\frac{1}{V}\int_{V}\mathbf{f}_{body}\mathrm{d}V.This covers the left-hand side of Newton’s second law, i.e.
\mathbf{F}=m\mathbf{a},now we need to work on the right-hand side. As we saw previously, we want to integrate over the entire system’s volume to account for each change in mass and acceleration, after all, that is exactly what a conservation law is there for, i.e. account for all contributions over a specific volume and track any changes. Therefore, we can already see that the mass will be divided by a volume which will provide us with units of kg/m^3, which, again, is the density. Therefore, we can already write
\int_V\mathbf{F}\mathrm{d}V=\int_V\rho\mathbf{a}\mathrm{d}V.Taking a units-based approach again, we can see that the acceleration a has units of m/s^2 and we can replace this term with the time derivative of the velocity, which is easier for us to measure. This can then be written as
\int_V\mathbf{F}\mathrm{d}V=\int_{V}\rho\frac{\mathrm{d}\mathbf{u}}{\mathrm{d} t}\mathrm{d}V.At this stage, we simply insert the definitions for the forces we developed above and apply the Reynolds Transport Theorem derived previously to the right-hand side, where we use \Phi=\rho\mathbf{u}. We also switch the equation around so that the time derivative is on the left-hand side, while all force terms are collected on the right-hand side. This results in
\int_{V}\frac{\partial \rho\mathbf{u}}{\partial t}\mathrm{d}V + \int_{V}\nabla\cdot (\rho\mathbf{u})\mathrm{d}V =
\int_{V}\nabla\cdot\sigma\mathrm{d}V + \int_{V}\rho\mathbf{g}\mathrm{d}V + \frac{1}{V}\int_{V}\mathbf{f}_{body}\mathrm{d}V.Carrying out the integration over the volume, we realise that the volume appears in each term and thus we can divide by it to arrive at the following, general momentum conservation equation (which is also known as the Cauchy momentum equation)
\frac{\partial \rho\mathbf{u}}{\partial t} + \rho(\mathbf{u}\cdot\nabla)\mathbf{u} = \nabla\cdot\sigma+\rho\mathbf{g}+\frac{1}{V}\mathbf{f}_{body}.The crucial aspect in deriving the momentum conservation laws for fluids is how we treat the internal forces, i.e. how we define \sigma. Without loss of generality, we can decompose it as alluded to above as a combination of the pressure and shear stresses as
\sigma = -p\mathbf{I} + \tau.The stress tensor \sigma has 9 components for a three-dimensional flow. However, the pressure forces can only act normal to surfaces (not tangential to it) and so we multiply it by the identity matrix \mathbf{I} to ensure its only acting normal to the x, y, and z direction and not in planes, which would have xy, xz, and yz components. The pressure is also a scalar value, acting in each direction with the same magnitude and the pressure itself is the same in each direction. The question then becomes what the shear stresses should look like. We have knowledge based on centuries of research on fluid mechanics, but that topic was not always that well understood.
Returning to our example of two blocks shearing against one another which produces a friction or shear force, if we were to neglect friction, that is to say, we are investigating a frictionless fluid, we could rewrite our stress tensor as
\sigma = -p\mathbf{I}.This is probably the simplest form of the Navier-Stokes equations we can come up with and this was used extensively throughout history and, when inserted into the Cauchy momentum equation, has become known as the Euler equation. The Euler equation was long accepted in the 19th century but the role of friction was not properly understood. Louis Marie Henri Navier, Jean Claude Barre de Saint-Venant, and George Gabriel Stokes all presented different approaches to extend the Euler equation to incorporate friction between fluid layers. In his 1822 paper titled “Memoire sur les lois du mouvement des fluides”, Navier presented his hypothesis about diffusion to the Academie des Sciences in Paris which was published five years later. His work was based on molecular considerations and interestingly, Anderson (J. D. Anderson, “Brief History of the Early Development of Theoretical and Experimental Fluid Dynamics,” 2010., p.12) highlights that:
“Navier’s equations were of the correct form, his theoretical reasoning was greatly flawed, and it is almost a fluke that he obtained the correct terms“.
It was not until the year 1843, seven years after the death of Navier, that Saint-Venant re-derived Navier’s hypothesis from a macroscopic point of view, completely ignoring the molecular considerations given by Navier. Around the same time, in 1845, Stokes published his hypothesis on diffusion in his paper titled “On the Theories of the Internal Friction of Fluids in Motion, and of the Equilibrium and Motion of Elastic Solids”. Stokes was unaware of the work of Navier or Saint-Venant. Despite the more rigorous derivation given by Saint-Venant, the equations became known as the Navier-Stokes equations.
It is worth noting here that if we inserted the definition for the stresses \sigma = -p\mathbf{I} + \tau. into our Cauchy momentum equation, we obtain
\frac{\partial \rho\mathbf{u}}{\partial t} + \rho(\mathbf{u}\cdot\nabla)\mathbf{u} = \nabla\cdot(-p\mathbf{I} + \tau)+\rho\mathbf{g}+\frac{1}{V}\mathbf{f}_{body}.The above equation is exact and there is no approximation done to it. The only issue arises, then, in determining the shear stresses \tau. If we take the same approach as above and set it to \tau=0, we obtain the inviscid Euler equations and these are still exact, so long as there are no viscous forces involved (which, in itself is a rather questionable assumption but for certain types of flows, e.g. compressible flows, this approximation provides useful results across a range of applications). What the work of Navier, Saint-Venant, and Stokes boils down to then, is how we impose viscous forces onto our momentum conservation law, that is to say, what form \tau should take in our the stress tensor \sigma = -p\mathbf{I} + \tau.
It is worth noting here that the assumption that Navier, Saint-Venant, and Stokes made for \tau results in a constitutive equation. This is an equation that can either be derived from experimental observations or from first principle, i.e. deriving it based on some fundamental physical laws (Newton’s laws, thermodynamics, conservation laws, etc.). The approximation these three put forward probably falls somewhere in between.
To understand how they came up with their constitutive equation, let’s assume that measuring the shear stress \tau would be easy for us to do in a lab, but if we wanted to find a relation for \tau, then we need to have some physical quantity to compare it against and put in relation to (this is exactly what a constitutive equation is providing, i.e. replacing an unknown quantity by something else we can measure). For example, if you wanted to measure the velocity of falling objects, how would you do that? Probably any mechanism you can come up with would involve some form of measuring the time it took for the object to fall a certain distance (video, laser, or sonar measurements). If we were then to plot the distances fallen for different time intervals, we would realise that there is some quadratic relationship between the time it takes for the object to fall a certain distance. And, from high-school physics, we can intuitively verify this as we know that a distance can be measured as x=0.5ut^2. How does that help us with our shear stress? Well, let’s try to think of another variable we can measure and then plot that variable against the shear stress itself. Trying to come up with a suitable variable is not as straightforward if you think about it. If we were to record the time or space at which point we measured the shear stresses, that would unlikely reveal any meaningful correlation. So, clever people before us have come up with hypotheses and tested these, among which Newton proposed that the shear stresses are related to the velocity gradient.
Before we look at a plot of shear stresses against the velocity gradient, we need to introduce a new physical property; the shear rate. The shear rate is given as \dot{\gamma}=u/h, where u is the relative velocity between two fluid layers in motion (typically taken between two solid plates in relative motion, i.e. the Couette flow) and h the distance between these two plates. Now if we let the distance between these two plates go to zero, noting that the relative velocity can be expressed as the difference between the fluid velocity at the upper and lower layer, then we can write
\dot{\gamma}=\frac{u}{h}=\lim_{h\rightarrow 0}\frac{u_{upper}-u_{lower}}{h}=\frac{\mathrm{d}u}{\mathrm{d}h}.If we replace the height h with the more general variable x, then we can see that the shear rate is nothing else than the velocity gradient in disguise. Thus, if we are able to measure the shear rate and the shear stresses, then we can plot their relationship and determine a mathematical model from that. As I mentioned above, countless people have done experiments on exactly that, and the below figure is reproduced from a study by B. C. Sahoo et al. 2009, where this relationship was measured in a lab for a nanofluid.
Inspecting the graph shows a linear relationship between the \tau and the shear rate (or velocity gradient) \dot{\gamma}. Thus, we can use basic linear equation fitting tools that we have probably picked up in high school as well, which follow the simple form
y=mx+b.
Here, we have y=\tau, x=\dot{\gamma}, and b=0, as all curves go through the origin at zero. We haven’t discussed the proportionality factor m here, but inspecting the above plot further, we can see that the only difference between these plots is the dynamic viscosity \mu, which is given for each curve in the graph. Thus, with m=\mu, we can rewrite the generic equation given above as
\tau=\mu\dot{\gamma}=\mu\frac{\mathrm{d}u}{\mathrm{d}y}The choice of the space coordinate along which we measure the velocity gradient here is arbitrary, but you may recognise the above form as the wall shear stresses. This concept plays an important role in modelling turbulence and is an omnipresent quantity in the study of fluid dynamics.
Now the above definition works well for cases where we have changes in velocity u along a flat plate, where the wall-normal direction is measured in y. This is typically the case in the turbulence literature, but in order to derive a general view of the Navier-Stokes equations and the stresses involved, we need to generalise this velocity gradient.
If we look back at our stress tensor, i.e. \sigma = -p\mathbf{I} + \tau, then we know that \sigma itself is a tensor, thus \tau needs to be a tensor as well and so we have to generalise the velocity gradient in our equation above for the shear stresses. A first assumption would be to simply write our velocity gradient tensor, i.e.
\nabla \mathbf{u}=\begin{bmatrix}\frac{\partial u}{\partial x} & \frac{\partial v}{\partial x} & \frac{\partial w}{\partial x} \\ \frac{\partial u}{\partial y} & \frac{\partial v}{\partial y} & \frac{\partial z}{\partial y} \\ \frac{\partial u}{\partial z} & \frac{\partial v}{\partial z} & \frac{\partial w}{\partial z}\end{bmatrix},and then use that definition in our shear stresses to generalise that expression to a tensor, i.e. \tau=\mu\nabla u. But if we think about this, the velocity gradient tensor does not represent stresses but just the derivatives of velocities. Thus, to understand what Navier, Saint-Venant, and Stokes did, we need to first decompose the velocity gradient tensor into a symmetric and anti-symmetric part as
\nabla \mathbf{u}=\frac{1}{2}\left[\nabla \mathbf{u}+(\nabla \mathbf{u})^T\right]+\frac{1}{2}\left[\nabla \mathbf{u}-(\nabla \mathbf{u})^T\right].The above decomposition can be done for any matrix or second-order tensor which was shown by the mathematician Toeplitz. To verify that this is correct, we can simply write out the tensors. introducing the shorthand notation \frac{\partial u}{\partial x}=u_x, we have
\nabla \mathbf{u}=\frac{1}{2}\begin{bmatrix}u_x & u_y & u_z \\ v_x & v_y & v_z \\ w_x & w_y & w_z\end{bmatrix}+\frac{1}{2}\begin{bmatrix}u_x & u_y & u_z \\ v_x & v_y & v_z \\ w_x & w_y & w_z\end{bmatrix}^T+\frac{1}{2}\begin{bmatrix}u_x & u_y & u_z \\ v_x & v_y & v_z \\ w_x & w_y & w_z\end{bmatrix}-\frac{1}{2}\begin{bmatrix}u_x & u_y & u_z \\ v_x & v_y & v_z \\ w_x & w_y & w_z\end{bmatrix}^T.We can see that the second and fourth tensor cancel out and we are left with
\nabla \mathbf{u}=\frac{1}{2}\begin{bmatrix}u_x & u_y & u_z \\ v_x & v_y & v_z \\ w_x & w_y & w_z\end{bmatrix}+\frac{1}{2}\begin{bmatrix}u_x & u_y & u_z \\ v_x & v_y & v_z \\ w_x & w_y & w_z\end{bmatrix}=\begin{bmatrix}u_x & u_y & u_z \\ v_x & v_y & v_z \\ w_x & w_y & w_z\end{bmatrix},which is the original velocity gradient tensor that we started with, confirming Toeplitz’s decomposition. Let’s investigate the anti-symmetric part (which is the third and forth tensor above), this can be written as
\frac{1}{2}\begin{bmatrix}u_x & u_y & u_z \\ v_x & v_y & v_z \\ w_x & w_y & w_z\end{bmatrix}-\frac{1}{2}\begin{bmatrix}u_x & v_x & w_x \\ u_y & v_y & w_y \\ u_z & v_z & w_z\end{bmatrix}^T = \frac{1}{2}\begin{bmatrix}u_x & u_y & u_z \\ v_x & v_y & v_z \\ w_x & w_y & w_z\end{bmatrix} - \frac{1}{2}\begin{bmatrix}u_x & v_x & w_x \\ u_y & v_y & w_y \\ u_z & v_z & w_z\end{bmatrix} =
\frac{1}{2}\begin{bmatrix}0 & u_y - v_x & u_z - w_x \\ v_x - u_y & 0 & v_z - w_y \\ w_x - u_z & w_y - v_z & 0\end{bmatrix}.Compare that to the curl of the velocity field
\nabla\times\mathbf{u}=\begin{bmatrix} w_y-v_z \\ u_z-w_x \\ v_x-u_y \end{bmatrix},and we can see how the tensor’s off-diagonal terms contain the same rotating velocity components, hence the anti-symmetric part is responsible for the rotation of the fluid. Rotation itself, however, does not contribute to any viscous losses. If you think about it, if we had two in-plane forces, say in the x and y direction, taking the cross product of those forces in these directions will result in a force in the z direction and therefore be normal to the plane spanned by the two forces in the x and y direction. Normal forces do not contribute to viscous shearing, which, by definition, can only happen in an in-plane direction, and thus Navier, Saint-Venant, and Stokes assumed the anti-symmetric part not to play a role in the viscous shear stress tensor.
Thus, they dropped it from the definition of the shear stresses and the first part of their hypothesis can be summarised as follows:
\tau=2\mu S=\mu[\nabla\mathbf{u}+(\nabla\mathbf{u})^T],where the symmetric part of the velocity gradient tensor is now abbreviated with the letter S, also referred to as the strain-rate tensor and we have combined it with Newton’s shear stress equation. We have to multiply by a factor of two as we have dropped the rotational part (anti-symmetric contribution). Since both parts contribute and equal half to the overall velocity gradient tensor, dropping one part means we have to multiply the other by a factor of two so as to not loose velocity gradient information.
In his 1845 paper On the theories of the internal friction of fluids in motion, George Gabriel Stokes derived his hypothesis of how friction in a fluid comes about and used molecular considerations to arrive at his final equations. He assumed that the velocity can be decomposed as
Total Velocity = Velocity due to Translation + Velocity due to Rotation + Velocity due to Dilatation
The first two terms are the ones we have seen above which were due to Toeplitz’s matrix decomposition into a symmetric and anti-symmetric part (note that the above is for the overall velocity, not shear stresses, hence we do account for the rotation here again). The third term, velocity due to dilatation requires a bit further explanation. Stokes used a small control volume that is moving at some point P in the flow which is moving with the mean velocities u, v and w. However, within that small control volume, if taken sufficiently small to resolve the molecular level, those molecules would have random velocity components. These velocity components would not have any preferred direction after the mean flow is removed (i.e. they are isotropically distributed) and so the collision of the molecules with the control volume would cause the control volume to expand (dilatation). This is shown below schematically.
He then postulated that this dilatation term should be
\delta=\nabla\cdot\mathbf{u}=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}For an incompressible flow, we have seen that the dilatation term is zero but for a general, compressible flow, it is not necessarily zero. This definition of the dilatation makes sense, say we have a 3D control volume and there is a change in the velocity u in the x direction, then either the second or third term above have to change equally in order for there not to be any dilatation. In other words, if the molecules in our control volume show a velocity gradient in one direction, that is to say, they are flowing in a preferred direction, our control volume needs to expand in that same direction in order to still contain all molecules. If the volume is changing, we have dilatation and thus this hypothesis makes intuitive sense.
Thus, Stokes retained the velocity due to translation but also added his third term due to dilatation to the definition of the stress tensor and arrived at
\tau = \mu[\nabla\mathbf{u}+(\nabla\mathbf{u})^T]+\lambda(\nabla\cdot\mathbf{u}).In order to understand the additional new variable \lambda, which is referred to as the bulk viscosity, Stokes needed to make some further assumptions. He assumed that the dilatation is uniform, i.e. the control volume is expanding or contracting in each direction uniformly. From a molecular point, this makes sense, as molecules have no preferred direction of movement once the bulk (mean) flow is removed from them. They show a stochastic, random motion that is equal in all directions and thus we could also write
\delta = \frac{1}{3}\nabla\cdot\mathbf{u}.This equation now assumes that all individual derivates \frac{\partial u}{\partial x}, \frac{\partial v}{\partial y} and \frac{\partial w}{\partial z} are equal and of the same magnitude and thus we need to divide by a factor of 3 in order to only take one of these directions into account (otherwise our dilatation would be 3 times as big as we were to account for it three times).
Next, Stokes also assumed that we can further write Newton’s shear stress equation as
\tau=2\mu[S-\delta],
Remember that the shear stresses are given in N/m^2, hence a positive dilatation would increase the control volume containing the molecules and thus we would have a larger surface area. If we have a larger surface area but the same shear stresses, a positive dilatation would then result in reduced shear stresses, hence we have to subtract the dilatation above from the strain rate tensor rather than adding the dilatation.
If we now insert Stokes’ hypothesis into our stress tensor, then we arrive at the final form of
\tau = \mu[\nabla\mathbf{u}+(\nabla\mathbf{u})^T]-\frac{2}{3}\mu(\nabla\cdot\mathbf{u}).At this point I should stress again that this is a constitutive equation, something these three people have come up with and we have found to work well for a range of applications, but that does not mean that they have provided us with a closure for the exact Cauchy momentum equation. Just because their mathematics looks elegant and we have found these equations to work well for a wide range of applications, we still have not definitely proved that these equations are indeed correct. In fact, there is still a lively debate as to whether the -\frac{2}{3}\mu term is correct or not. The uniqueness and existence of the Navier-Stokes equations remains an open problem and is one of the 7 unsolved problems listed by the Clay Institute, for which there is a 1 000 000 $ award to be claimed by someone able to provide proof of these equations. You can read up on this on the Clay Institute website (and post your solution in the comment below, sharing is caring, right? …).
Taking the derivative of the above equation, the following result is obtained
\nabla\cdot\tau = \mu[\nabla\cdot\nabla\mathbf{u} + \nabla(\nabla\cdot\mathbf{u})] - \frac{2}{3}\mu\nabla(\nabla\cdot\mathbf{u}).Inserting this definition of the viscous stress tensor into \sigma = -p\mathbf{I} + \tau and then into the Cauchy momentum equation, and dividing all terms by the density, the following equation is obtained
\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \nu\nabla(\nabla\cdot\mathbf{u}) - \frac{2}{3}\nu\nabla(\nabla\cdot\mathbf{u}) + \mathbf{g} +\frac{1}{\rho V}\mathbf{f}_{body},where \nu=\mu/\rho is the kinematic viscosity. Defining Q=\nu\nabla(\nabla\cdot\mathbf{u}) - \frac{2}{3}\nu\nabla(\nabla\cdot\mathbf{u}) + \mathbf{g} +\frac{1}{\rho V}\mathbf{f}_{body} as a source term, we can write the above equation as
\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + Q,This is the conservation of momentum equation according to Navier, Saint-Venant, and Stokes, or, as we call it, the Navier-Stokes equation. Originally the Navier-Stokes equation was simply that, i.e. the conservation of momentum, but these days we also consider the conservation of mass and energy to be part of the set of Navier-Stokes equations, although Navier and Stokes had nothing to do with these conservation laws.
At this point it is worthwhile to point out that the incompressible Navier-Stokes equations are obtained by realising that we got \nabla\cdot\mathbf{u}=0 from the continuity equation, and so if we insert that into our generalised Navier-Stokes momentum equation, we obtain
\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \mathbf{g} +\frac{1}{\rho V}\mathbf{f}_{body}.Most of the time this equation is simplified further by dropping body force terms and the influence of gravity, and we obtain the following simplified equation
\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u}.Summary
To summarise what we have derived above, the following equations were obtained.
The Reynolds Transport Theorem
\frac{\mathrm{d}\Phi(t)}{\mathrm{d}t}=\frac{\partial \phi}{\partial t}+\nabla\cdot\phi\mathbf{u}.Conservation of Mass (compressible flow)
\frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\mathbf{u})=0.Conservation of Mass (incompressible flow)
\nabla\cdot\mathbf{u}=0Conservation of Momentum (compressible flow)
\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \nu\nabla(\nabla\cdot\mathbf{u}) - \frac{2}{3}\nu\nabla(\nabla\cdot\mathbf{u}) + \mathbf{g} +\frac{1}{\rho V}\mathbf{f}_{body},Conservation of Momentum (incompressible flow)
\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \mathbf{g} +\frac{1}{\rho V}\mathbf{f}_{body}.Before You go
I started this website to build and bring together a community of like-minded people in the field of CFD. If you liked the article and would like to be part of that, then consider signing up to my newsletter. I won’t spam you but would love to see you on board. If you are ready, let me know by signing up!