Derivation of Navier-Stokes Equation in cylindrical coordinates

The Navier-Stokes equations are derived in cylindrical coordinates. Since the strain rate tensor $\mathbf{D}=(\mathbf{L}+{\mathbf{L}}^T)/2$, if $\mathbf{V}=v_r{\mathbf{e}}_{\mathbf{r}}+v_{\theta }{\mathbf{v}}_{\theta }+v_z{\mathbf{e}}_{\mathbf{z}}$ (where $L=\mathrm{\nabla }V$ is the gradient tensor of velocity vector), the strain rate and stress tensors are

$$\begin{array}{rl}& \mathbf{D}=\left[\begin{array}{ccc}{\displaystyle \frac{\partial v_r}{\partial r}}& {\displaystyle \frac{1}{2}}({\displaystyle \frac{\partial v_r}{r\partial \theta }}-{\displaystyle \frac{v_{\theta }}{r}}+{\displaystyle \frac{\partial v_{\theta }}{\partial r}})& {\displaystyle \frac{1}{2}}({\displaystyle \frac{\partial v_r}{\partial z}}+{\displaystyle \frac{\partial v_z}{\partial r}})\\ {\displaystyle \frac{1}{2}}({\displaystyle \frac{\partial v_r}{r\partial \theta }}-{\displaystyle \frac{v_{\theta }}{r}}+{\displaystyle \frac{\partial v_{\theta }}{\partial r}})& {\displaystyle \frac{\partial v_{\theta }}{r\partial \theta }}+{\displaystyle \frac{v_r}{r}}& {\displaystyle \frac{1}{2}}({\displaystyle \frac{\partial v_{\theta }}{\partial z}}+{\displaystyle \frac{\partial v_z}{r\partial \theta }})\\ {\displaystyle \frac{1}{2}}({\displaystyle \frac{\partial v_r}{\partial z}}+{\displaystyle \frac{\partial v_z}{\partial r}})& {\displaystyle \frac{1}{2}}({\displaystyle \frac{\partial v_{\theta }}{\partial z}}+{\displaystyle \frac{\partial v_z}{r\partial \theta }})& {\displaystyle \frac{\partial v_z}{\partial z}}\end{array}\right]\\ & \mathit{\sigma }=\left[\begin{array}{ccc}{\sigma }_{rr}& {\tau }_{r\theta }& {\tau }_{rz}\\ {\tau }_{\theta r}& {\sigma }_{\theta \theta }& {\tau }_{\theta z}\\ {\tau }_{zr}& {\tau }_{z\theta }& {\sigma }_{zz}\end{array}\right]\end{array}$$

By substituting these matices into equation (1) in the post of 'Navier-Stokes in cartesian coordinates' and comparing each component, the relationship between the stresses and the velocity gradients is as follows.

$$\begin{array}{rl}& {\sigma }_{rr}=-p-\frac{2}{3}\mu \mathrm{\nabla }\cdot \mathbf{V}+2\mu \frac{\partial v_r}{\partial r}\\ & {\sigma }_{\theta \theta }=-p-\frac{2}{3}\mu \mathrm{\nabla }\cdot \mathbf{V}+2\mu (\frac{\partial v_{\theta }}{r\partial \theta }+\frac{v_r}{r})\\ & {\sigma }_{zz}=-p-\frac{2}{3}\mu \mathrm{\nabla }\cdot \mathbf{V}+2\mu \frac{\partial v_z}{\partial z}\\ & {\tau }_{r\theta }={\tau }_{\theta r}=\mu \{\frac{r\partial }{\partial r}\left(\frac{v_{\theta }}{r}\right)+\frac{\partial v_r}{r\partial \theta }\}\\ & {\tau }_{rz}={\tau }_{zr}=\mu (\frac{\partial v_r}{\partial z}+\frac{\partial v_z}{\partial r})\\ & {\tau }_{\theta z}={\tau }_{z\theta }=\mu (\frac{\partial v_{\theta }}{\partial z}+\frac{\partial v_z}{r\partial \theta })\end{array}$$

And, the equation of motion is given as follows by substituting the Eulerian expression for acceleration, $D\mathbf{V}/Dt=\mathbf{a}$, and the Laplacian of the velocity vector, ${\mathrm{\nabla }}^2\mathbf{V}$, into equation (2) in the cylindrical coordinate system.

$$\begin{array}{rl}& \frac{\partial v_r}{\partial t}+v_r\frac{\partial v_r}{\partial r}+v_{\theta }\frac{\partial v_r}{r\partial \theta }+v_z\frac{\partial v_r}{\partial z}-\frac{v_{\theta }^2}{r}\\ & =-\frac{\partial p}{\rho \partial r}+\frac{\nu \partial }{3\partial r}(\mathrm{\nabla }\cdot \mathbf{V})+\nu \{\frac{\partial }{r\partial r}\left(r\frac{\partial v_r}{\partial r}\right)+\frac{{\partial }^2{v}_r}{r^2\partial {\theta }^2}+\frac{{\partial }^2{v}_r}{\partial z^2}-\frac{v_r}{r^2}-\frac{2\partial v_{\theta }}{r^2\partial \theta }\}+f_r\\ & \frac{\partial v_{\theta }}{\partial t}+v_r\frac{\partial v_{\theta }}{\partial r}+v_{\theta }\frac{\partial v_{\theta }}{r\partial \theta }+v_z\frac{\partial v_{\theta }}{\partial z}+\frac{v_r{v}_{\theta }}{r}\\ & =-\frac{\partial p}{\rho r\partial \theta }+\frac{\nu \partial }{3\partial r}(\mathrm{\nabla }\cdot \mathbf{V})+\nu \{\frac{\partial }{r\partial r}\left(r\frac{\partial v_{\theta }}{\partial r}\right)+\frac{{\partial }^2{v}_{\theta }}{r^2\partial {\theta }^2}+\frac{{\partial }^2{v}_{\theta }}{\partial z^2}-\frac{v_{\theta }}{r^2}+\frac{2\partial v_r}{r^2\partial \theta }\}+f_{\theta }\\ & \frac{\partial v_z}{\partial t}+v_r\frac{\partial v_z}{\partial r}+v_{\theta }\frac{\partial v_z}{r\partial \theta }+v_z\frac{\partial v_z}{r\partial z}\\ & =-\frac{\partial p}{\rho \partial z}+\frac{\nu \partial }{3\partial z}(\mathrm{\nabla }\cdot \mathbf{V})+\nu \{\frac{\partial }{r\partial r}\left(r\frac{\partial v_z}{\partial r}\right)+\frac{{\partial }^2{v}_z}{r^2\partial z^2}+\frac{{\partial }^2{v}_z}{\partial z^2}\}+f_z\end{array}$$

[Example] Find the equations of motion for incompressible, laminar, steady flow in a circular tube, assuming that the flow has only an axial(z-direction) velocity component.

[Solution] The equations of motion for floq within a tube need to be expressed in cylindrical coordinates. There are two methods for obtaining the equations of motion: one is to set up a small element in the flow field and directly apply Newton's laws of motion to it, and the other is to simplify the general equations of motion, i.e. the Navier-Stokes equations. In this example, both of these methods will be explained.

(1) Since the given flow is axisymmetric, we set up an annular element like the figure in the flow field. If the fluid is not accelerated, the relationship of
$$d\mathbf{F}=dm\mathbf{a}=0\cdots (\mathrm{a})$$

will be established, therefore the shear force acting on the element will be balanced by the pressure difference from the pressure gradient, $dp/dz$. These forces are shown in the figure. If we first consider the equilibrium of the axial force, it it as follows.

$$\begin{array}{rl}{F}_z& =p2\pi (r+\frac{\mathrm{\Delta }r}{2})\mathrm{\Delta }r-(p+\frac{\partial p}{\partial z}\mathrm{\Delta }z)2\pi (r+\frac{\mathrm{\Delta }r}{2})\mathrm{\Delta }r-{\tau }_{rz}2\pi r\mathrm{\Delta }z\\ & +({\tau }_{rz}+\frac{\partial {\tau }_{rz}}{\partial r}\mathrm{\Delta }r)2\pi (r+\mathrm{\Delta }r)\mathrm{\Delta }z\\ & =-\frac{\partial p}{\partial z}\mathrm{\Delta }z2\pi r\mathrm{\Delta }r-\frac{\partial p}{\partial z}\mathrm{\Delta }z\pi \mathrm{\Delta }{r}^2+{\tau }_{rz}2\pi \mathrm{\Delta }r\mathrm{\Delta }z+\frac{\partial {\tau }_{rz}}{\partial r}\mathrm{\Delta }r2\pi r\mathrm{\Delta }z=0\cdots (\mathrm{b})\end{array}$$

Figure   Annular infinitesimal volume

Where, when $\mathrm{\Delta }r$ and $\mathrm{\Delta }z$ have values close to zero($\mathrm{\Delta }{r}^2\approx 0$), the above equation is summarized as follow.

$$-\frac{\partial p}{\partial z}+\frac{{\tau }_{rz}}{r}+\frac{\partial {\tau }_{rz}}{\partial r}=0\cdots (\mathrm{c})$$

However, in the case of Newtonian fluid,

$${\tau }_{rz}=\mu \frac{\partial v_z}{\partial r}\cdots (\mathrm{d})$$

can be written, then substitute it into the equation (c), we can obtain

$$-\frac{\partial p}{\partial z}+\frac{\mu \partial v_z}{r\partial r}+\mu \frac{{\partial }^2{v}_z}{\partial r^2}=0\cdots (\mathrm{e})$$

The given flow is axisymmetric, hence no change in $\theta$ direction($\partial /\partial \theta =0$), from the equilibrium of force in radial direction(r direction),

$$\frac{\partial p}{\partial r}=0\cdots (\mathrm{f})$$

can be obtained, hence pressure is a function of only $z$. In addition, since $v_z$ is also a function of $r$ alone, the following equation of motion is obtained.

$$-\frac{dp}{dz}+\frac{\mu d}{rdr}\left(r\frac{d{v}_z}{dr}\right)=0\cdots (\mathrm{g})$$

(2) From the general continuity equation, for incompressible flow, we obtain the following equation.

$$\mathrm{\nabla }\cdot \mathrm{V}=\frac{\partial (r{v}_r)}{r\partial r}+\frac{\partial v_{\theta }}{r\partial \theta }+\frac{\partial v_z}{\partial z}=0\cdots (\mathrm{h})$$

However, since the given flow has only the velocity component in the $z$ direction, therefore

$$v_r=v_{\theta }=0\cdots (\mathrm{i})$$

Substituting this into the equation (h), we can get the relationship of

$$\frac{\partial v_z}{\partial z}=0\cdots (\mathrm{j})$$

And $\partial /\partial \theta =0$ in an axisymmetric flow. By substituting all of these conditions into the equations of motion in cylindrical coordinates, we get the equations below.

$$\begin{array}{rl}& \frac{\partial p}{\partial r}=0\cdots (\mathrm{k})\\ & \frac{\partial p}{\partial \theta }=0\cdots (\mathrm{l})\\ -& \frac{\partial p}{\partial z}+\frac{\mu \partial }{r\partial r}\left(r\frac{\partial v_z}{\partial r}\right)=0\cdots (\mathrm{m})\end{array}$$

Equation (k) and (l) show that the pressure, $p$ is a function of $z$ only, and from equation (j) and the axisymmetric condition, we get to know that $v_z$ is also a function of $r$ only. Therefore, we get

$$-\frac{dp}{dz}+\frac{\mu d}{rdr}\left(r\frac{d{v}_z}{dr}\right)=0\cdots (\mathrm{n})$$

This equation is equal to the equation (g) obtained previously.

As shown in this example, the equations of motion for a given flow field can be obtained more easily by using the general equations of motion, the Navier-Stokes equation.