Stresses in the plane orthogonal to \(r\) and \(z\) direction. This figure is taken from https://scribd.com/doc/111812009/Navier-Stokes-Derivation-in-Cylindrical-Coordiantes.
The Cauchy momentum equation is a vector partial differential equation put forth by Cauchy that describes the non-relativistic momentum transport in any continuum. It expresses Newton’s second law and in convective (or Lagrangian) form it can be written as: \[\label{Eq:cauchy} \frac{D\mathbf{u}}{Dt}=\frac{1}{\rho}\nabla\cdot {\boldsymbol{\mathbf{\sigma}}}+\mathbf{g}\] where \(\rho\) is the density at the point considered in the continuum (for which the continuity equation holds), \({\boldsymbol{\mathbf{\sigma}}}\) is the stress tensor, and \(\mathbf{g}\) contains all of the body forces per unit mass (often simply gravitational acceleration). In the Cauchy equation \(\mathbf{u}\) is the flow velocity vector field, which depends on time and space.
We want to write the terms of Eq. (1) in cylindrical coordinates. First of all, we write the flow velocity vector in cylindrical coordinates as: \[\mathbf{u}(r,\theta,z,t) = u_r(r,\theta,z,t)\mathbf{e}_r + u_{\theta}(r,\theta,z,t)\mathbf{e}_{\theta} + u_z(r,\theta,z,t)\mathbf{e}_z,\] where \(\left\{ {\boldsymbol{\mathbf{e}}}_r,{\boldsymbol{\mathbf{e}}}_{\theta},{\boldsymbol{\mathbf{e}}}_z \right\}\) is a right-handed triad of unit vectors.
Let consider first the material derivative on the left-hand side: \[\label{Eq:material_derivative} \frac{D\mathbf{u}}{Dt}=\frac{\partial {\mathbf{u}}}{\partial t} + \frac{\partial {\mathbf{u}}}{\partial r} \frac{d r}{dt} + \frac{\partial {\mathbf{u}}}{\partial \theta} \frac{d \theta}{dt} + \frac{\partial {\mathbf{u}}}{\partial z} \frac{d z}{dt}.\]
When differentiating the velocity vector in cylindrical coordinates, the unit vectors must also be differentiated, because they are not fixed. In particular, we will make use of the following relations: \[\label{eq:versor_deriv} \frac{\partial {\boldsymbol{\mathbf{e}}}_{r}}{\partial r} = \frac{\partial {\boldsymbol{\mathbf{e}}}_{\theta}}{\partial r} = \frac{\partial {\boldsymbol{\mathbf{e}}}_{z}}{\partial r} = 0, \quad \frac{\partial {\boldsymbol{\mathbf{e}}}_{r}}{\partial \theta} = {\boldsymbol{\mathbf{e}}}_{\theta}, \quad \frac{\partial {\boldsymbol{\mathbf{e}}}_{\theta}}{\partial \theta} = - {\boldsymbol{\mathbf{e}}}_{r}.\]
The first partial derivative on the right-hand side of Eq. (3) is computed at a fixed position in time, thus the unit vectors do not change in time and theur derivatives are identically zero. Then, we have \[\frac{\partial {\mathbf{u}}}{\partial t} = \frac{\partial u_r}{\partial t} \mathbf{e}_r + \frac{\partial u_{\theta}}{\partial t} \mathbf{e}_{\theta} + \frac{\partial u_z}{\partial t} \mathbf{e}_z.\]
By making use of the definition of velocities (\(u_r=dr/dt\), \(u_{\theta}/r=d\theta/dt\) and \(u_z=dz/dt\)), we can now evaluate the remaining terms of Eq. (3) as follows
\[\begin{array}{rl} \frac{dr}{dt}\frac{\partial \mathbf{u}}{\partial r} & = u_r \frac{\partial \mathbf{u}}{\partial r} \\ \\ & = u_r \left[ \frac{\partial u_r}{\partial r} \mathbf{e}_r + u_r \frac{\partial \mathbf{e}_r}{\partial r} + \frac{\partial u_{\theta}}{\partial r} \mathbf{e}_{\theta} + u_{\theta} \frac{\partial \mathbf{e}_{\theta}}{\partial r} + \frac{\partial u_z}{\partial r} \mathbf{e}_z + u_z \frac{\partial \mathbf{e}_z}{\partial r}\right] \\ \\ & = u_r \left[ \frac{\partial u_r}{\partial r} \mathbf{e}_r + \frac{\partial u_{\theta}}{\partial r} \mathbf{e}_{\theta} + \frac{\partial u_z}{\partial r} \mathbf{e}_z \right] \end{array}\]
\[\begin{array}{rl} \frac{d{\theta}}{dt}\frac{\partial \mathbf{u}}{\partial {\theta}} & = \frac{ u_{\theta} }{r} \frac{\partial \mathbf{u}}{\partial {\theta}} \\ \\ & = \frac{ u_{\theta} }{r} \left[ \frac{\partial u_r}{\partial {\theta}} \mathbf{e}_r + u_r \frac{\partial \mathbf{e}_r}{\partial {\theta}} + \frac{\partial u_{\theta}}{\partial {\theta}} \mathbf{e}_{\theta} + u_{\theta} \frac{\partial \mathbf{e}_{\theta}}{\partial {\theta}} + \frac{\partial u_z}{\partial {\theta}} \mathbf{e}_z + u_z \frac{\partial \mathbf{e}_z}{\partial {\theta}}\right] \\ \\ & = \frac{ u_{\theta} }{r} \left[ \frac{\partial u_r}{\partial {\theta}} \mathbf{e}_r + u_r \mathbf{e}_{\theta} + \frac{\partial u_{\theta}}{\partial {\theta}} \mathbf{e}_{\theta} + u_{\theta} (-\mathbf{e}_r) + \frac{\partial u_z}{\partial {\theta}} \mathbf{e}_z \right] \\ \\ & = \frac{ u_{\theta} }{r} \left[ ( \frac{\partial u_r}{\partial {\theta}} -u_{\theta} ) \mathbf{e}_r + ( u_r + \frac{\partial u_{\theta}}{\partial {\theta}} ) \mathbf{e}_{\theta} + \frac{\partial u_z}{\partial {\theta}} \mathbf{e}_z \right] \end{array}\]
\[\begin{array}{rl} \frac{dz}{dt}\frac{\partial \mathbf{u}}{\partial z} & = u_z \frac{\partial \mathbf{u}}{\partial z} \\ \\ & = u_z \left[ \frac{\partial u_z}{\partial z} \mathbf{e}_r + u_z \frac{\partial \mathbf{e}_r}{\partial z} + \frac{\partial u_{\theta}}{\partial z} \mathbf{e}_{\theta} + u_{\theta} \frac{\partial \mathbf{e}_{\theta}}{\partial z} + \frac{\partial u_z}{\partial z} \mathbf{e}_z + u_z \frac{\partial \mathbf{e}_z}{\partial z}\right] \\ \\ & = u_z \left[ \frac{\partial u_r}{\partial z} \mathbf{e}_r + \frac{\partial u_{\theta}}{\partial z} \mathbf{e}_{\theta} + \frac{\partial u_z}{\partial z} \mathbf{e}_z \right] \end{array}\]
Putting all the terms together and collecting the terms multiplying the same unit vector we can write the material derivative as:
\[\label{eq:material_expanded} \begin{array}{rl} \frac{D\mathbf{u}}{Dt} &= \left[ \frac{\partial u_r}{\partial t} + u_r \frac{\partial u_r}{\partial r} + \frac{ u_{\theta}} {r} \frac{\partial u_r}{\partial \theta} - \frac{u_{\theta}^2}{r} + u_z \frac{\partial u_r}{\partial z} \right] \mathbf{e}_r \\ \\ &+ \left[ \frac{\partial u_{\theta}}{\partial t} + u_r \frac{\partial u_{\theta}}{\partial r} + \frac{ u_{\theta}} {r} \frac{\partial u_{\theta}}{\partial \theta} + \frac{u_{\theta}u_r}{r} + u_z \frac{\partial u_{\theta}}{\partial z} \right] \mathbf{e}_{\theta} \\ \\ &+ \left[ \frac{\partial u_z}{\partial t} + u_r \frac{\partial u_z}{\partial r} + \frac{ u_{\theta}} {r} \frac{\partial u_z}{\partial \theta} + u_z \frac{\partial u_z}{\partial z} \right] \mathbf{e}_z. \end{array}\]
We want to evaluate here the term \(\nabla\cdot{\boldsymbol{\mathbf{\sigma}}}\) appearing in the Cauchy momentum equation in cylindrical coordinates. To this aim we compute the term for an infinitesimal volume as represented in Figure 1.
Stresses in the plane orthogonal to \(r\) and \(z\) direction. This figure is taken from https://scribd.com/doc/111812009/Navier-Stokes-Derivation-in-Cylindrical-Coordiantes.
By summing all the contributions in the \(r\) direction given by the stresses at the faces of the control volume we have:
\[\begin{array}{l} \left( \sigma_{rr} + \frac{\partial \sigma_{rr}}{\partial r} \delta r \right) (r+\delta r) \delta \theta \delta z - \sigma_{rr} r \delta \theta \delta z \\ \\ + \left( \sigma_{r\theta} + \frac{\partial \sigma_{r\theta}}{\partial \theta} \delta \theta \right) \delta r \delta z \cos \frac{\delta\theta}{2} - \sigma_{r\theta} \delta r \delta z \cos \frac{\delta\theta}{2} \\ \\ + \left( \sigma_{rz} + \frac{\partial \sigma_{rz}}{\partial z} \delta z \right) \left( r + \frac{\delta r}{2} \right) \delta r \delta \theta - \sigma_{rz} \left( r + \frac{\delta r}{2} \right) \delta r \delta \theta \\ \\ - \left( \sigma_{\theta\theta} + \frac{\partial \sigma_{\theta\theta}}{\partial \theta} \delta \theta \right) \delta r \delta z \sin \frac{\delta\theta}{2} - \sigma_{\theta\theta} \delta r \delta z \sin \frac{\delta\theta}{2} \end{array}\]
Canceling some terms, dividing by the volume \(r\delta r\delta\theta \delta z\), and letting \(\delta r,\delta\theta\) and \(\delta z\) go to zero, it yields:
\[\frac{\sigma_{rr}}{r} + \frac{\partial \sigma_{rr}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{r\theta}}{\partial \theta} + \frac{\partial \sigma_{rz}}{\partial z} - \frac{\sigma_{\theta\theta}}{r}.\]
Similarly, we find for the \(\theta\) direction: \[\frac{2\sigma_{r\theta}}{r} + \frac{\partial \sigma_{r\theta}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{\theta\theta}}{\partial \theta} + \frac{\partial \sigma_{\theta z}}{\partial z},\] and finally, for the \(z\) direction: \[\frac{\sigma_{rz}}{r} + \frac{\partial \sigma_{rz}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{\theta z}}{\partial \theta} + \frac{\partial \sigma_{z z}}{\partial z}.\] We can thus write the stress term in Eq. (1) in the following way: \[\label{Eq:stress_expanded} \begin{array}{rl} \nabla\cdot {\boldsymbol{\mathbf{\sigma}}} &= \left( \frac{\sigma_{rr}}{r} + \frac{\partial \sigma_{rr}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{r\theta}}{\partial \theta} + \frac{\partial \sigma_{rz}}{\partial z} - \frac{\sigma_{\theta\theta}}{r} \right) \mathbf{e}_r\\ \\ &+ \left( \frac{2\sigma_{r\theta}}{r} + \frac{\partial \sigma_{r\theta}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{\theta\theta}}{\partial \theta} + \frac{\partial \sigma_{\theta z}}{\partial z} \right) \mathbf{e}_{\theta}\\ \\ &+ \left( \frac{\sigma_{rz}}{r} + \frac{\partial \sigma_{rz}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{\theta z}}{\partial \theta} + \frac{\partial \sigma_{z z}}{\partial z} \right) \mathbf{e}_z. \end{array}\]
The previous expression could also be obtained writing the stress tensor \({\boldsymbol{\mathbf{\sigma}}}\) with the following dyadic representation: \[\begin{array}{rl} {\boldsymbol{\mathbf{\sigma}}} = & \sigma_{rr} {\boldsymbol{\mathbf{e}}}_r \otimes {\boldsymbol{\mathbf{e}}}_{r} + \sigma_{r\theta} {\boldsymbol{\mathbf{e}}}_r \otimes {\boldsymbol{\mathbf{e}}}_{\theta} + \sigma_{rz} {\boldsymbol{\mathbf{e}}}_r \otimes {\boldsymbol{\mathbf{e}}}_{z} \\ & + \sigma_{\theta r} {\boldsymbol{\mathbf{e}}}_{\theta} \otimes {\boldsymbol{\mathbf{e}}}_{r} + \sigma_{\theta\theta} {\boldsymbol{\mathbf{e}}}_{\theta} \otimes {\boldsymbol{\mathbf{e}}}_{\theta} + \sigma_{\theta z} {\boldsymbol{\mathbf{e}}}_{\theta} \otimes {\boldsymbol{\mathbf{e}}}_{z} \\ & + \sigma_{z r} {\boldsymbol{\mathbf{e}}}_{z} \otimes {\boldsymbol{\mathbf{e}}}_{r} + \sigma_{z\theta} {\boldsymbol{\mathbf{e}}}_{z} \otimes {\boldsymbol{\mathbf{e}}}_{\theta} + \sigma_{zz} {\boldsymbol{\mathbf{e}}}_{z} \otimes {\boldsymbol{\mathbf{e}}}_{z} \end{array}\]
The divergence of \({\boldsymbol{\mathbf{\sigma}}}\) is a vector, which can be represented as the operator1 \[\left( {\boldsymbol{\mathbf{e}}}_{r} \frac{\partial}{\partial r} + {\boldsymbol{\mathbf{e}}}_{\theta} \frac{1}{r}\frac{\partial}{\partial \theta} + {\boldsymbol{\mathbf{e}}}_{z} \frac{\partial}{\partial z}\right)\cdot \left( \begin{array}{c} \sigma_{rr} {\boldsymbol{\mathbf{e}}}_r \otimes {\boldsymbol{\mathbf{e}}}_{r} + \sigma_{r\theta} {\boldsymbol{\mathbf{e}}}_r \otimes {\boldsymbol{\mathbf{e}}}_{\theta} + \sigma_{rz} {\boldsymbol{\mathbf{e}}}_r \otimes {\boldsymbol{\mathbf{e}}}_{z} \\ + \sigma_{\theta r} {\boldsymbol{\mathbf{e}}}_{\theta} \otimes {\boldsymbol{\mathbf{e}}}_{r} + \sigma_{\theta\theta} {\boldsymbol{\mathbf{e}}}_{\theta} \otimes {\boldsymbol{\mathbf{e}}}_{\theta} + \sigma_{\theta z} {\boldsymbol{\mathbf{e}}}_{\theta} \otimes {\boldsymbol{\mathbf{e}}}_{z} \\ + \sigma_{z r} {\boldsymbol{\mathbf{e}}}_{z} \otimes {\boldsymbol{\mathbf{e}}}_{r} + \sigma_{z\theta} {\boldsymbol{\mathbf{e}}}_{z} \otimes {\boldsymbol{\mathbf{e}}}_{\theta} + \sigma_{zz} {\boldsymbol{\mathbf{e}}}_{z} \otimes {\boldsymbol{\mathbf{e}}}_{z} \end{array} \right)\]
Evaluating the components of the divergence of a tensor is an extremely tedious operation because, in addition to the components themselves, each of the basis vectors in the dyadic representation of \({\boldsymbol{\mathbf{\sigma}}}\) must be differentiated (see Eq. (4). For this reason I do not present the full derivation but only the evaluation of terms of the previous expression that contribute to the \(z\)-component of the term \(\nabla\cdot {\boldsymbol{\mathbf{\sigma}}}\):
\[\begin{aligned} &\left( {\boldsymbol{\mathbf{e}}}_{r} \frac{\partial}{\partial r} \right)\cdot \left( \sigma_{rz} {\boldsymbol{\mathbf{e}}}_r \otimes {\boldsymbol{\mathbf{e}}}_{z} \right) + \left( {\boldsymbol{\mathbf{e}}}_{\theta} \frac{1}{r} \frac{\partial}{\partial \theta} \right)\cdot \left( \sigma_{r z} {\boldsymbol{\mathbf{e}}}_r \otimes {\boldsymbol{\mathbf{e}}}_{z} + \sigma_{\theta z} {\boldsymbol{\mathbf{e}}}_{\theta} \otimes {\boldsymbol{\mathbf{e}}}_{z} + \right) \\ \\ &+ \left( {\boldsymbol{\mathbf{e}}}_{z} \frac{\partial}{\partial z} \right)\cdot \left( \sigma_{zz} {\boldsymbol{\mathbf{e}}}_r \otimes {\boldsymbol{\mathbf{e}}}_{z} \right)= \\ \\ & \frac{\partial \sigma_{rz}}{\partial r} ({\boldsymbol{\mathbf{e}}}_r \cdot {\boldsymbol{\mathbf{e}}}_r) {\boldsymbol{\mathbf{e}}}_z + \frac{\sigma_{rz}}{r} \left( \frac{\partial {\boldsymbol{\mathbf{e}}}_r}{\partial \theta} \cdot {\boldsymbol{\mathbf{e}}}_\theta \right){\boldsymbol{\mathbf{e}}}_z + \\ \\ & \frac{1}{r} \frac{\partial \sigma_{\theta z}}{\partial \theta} ({\boldsymbol{\mathbf{e}}}_{\theta} \cdot {\boldsymbol{\mathbf{e}}}_{\theta}) {\boldsymbol{\mathbf{e}}}_z + \frac{\partial \sigma_{zz}}{\partial z} ({\boldsymbol{\mathbf{e}}}_z \cdot {\boldsymbol{\mathbf{e}}}_z) {\boldsymbol{\mathbf{e}}}_z = \\ \\ & \left( \frac{\partial \sigma_{rz}}{\partial r} + \frac{\sigma_{rz}}{r} + \frac{1}{r} \frac{\partial \sigma_{\theta z}}{\partial \theta} + \frac{\partial \sigma_{z z}}{\partial z} \right) {\boldsymbol{\mathbf{e}}}_z\end{aligned}\]
The final result coincides with that presented in Eq. ([Eq:stress_expanded]).
The effect of stress in the continuum flow can be represented by the \(\nabla p\) and \(\nabla\cdot {\boldsymbol{\mathbf{\tau}}}\) terms; these are gradients of surface forces, analogous to stresses in a solid. Here \(\nabla p\) is called the pressure gradient and arises from the isotropic part of the Cauchy stress tensor, which has order two. This part is given by normal stresses that turn up in almost all situations, dynamic or not. The anisotropic part of the stress tensor gives rise to \(\nabla\cdot {\boldsymbol{\mathbf{\tau}}}\), which conventionally describes viscous forces; for incompressible flow, this is only a shear effect. Thus, \({\boldsymbol{\mathbf{\tau}}}\) is the deviatoric stress tensor, and the stress tensor is equal to: \[\label{eq:stress_tensor} {\boldsymbol{\mathbf{\sigma}}}=-p\mathbf{Id}+{\boldsymbol{\mathbf{\tau}}},\] where \(\mathbf{Id}\) is the identity matrix in the space considered and \({\boldsymbol{\mathbf{\tau}}}\) the shear tensor.
For an incompressible fluid with constant viscosity, the shear tensor can be written as \[\label{eq:stress_strain} {\boldsymbol{\mathbf{\tau}}} = \mu \left( \nabla \mathbf{u} + \nabla \mathbf{u}^T\right).\]
The gradient of the vector \({\boldsymbol{\mathbf{u}}}\) is a tensor, which can be represented as a dyadic product of the vector with the gradient operator as \[\begin{aligned} {\boldsymbol{\mathbf{u}}}\otimes\nabla = \left( u_r {\boldsymbol{\mathbf{e}}}_{r} + u_{\theta} {\boldsymbol{\mathbf{e}}}_{\theta} + u_z {\boldsymbol{\mathbf{e}}}_{z} \right) \otimes \left( {\boldsymbol{\mathbf{e}}}_{r} \frac{\partial}{\partial r} + {\boldsymbol{\mathbf{e}}}_{\theta} \frac{1}{r}\frac{\partial}{\partial \theta} + {\boldsymbol{\mathbf{e}}}_{z} \frac{\partial}{\partial z}\right).\end{aligned}\] The dyadic product can be expanded, but when evaluating the derivatives it is important to recall that the basis vectors are functions of the coordinate \(\theta\) and consequently their derivatives may not vanish. For example \[\begin{aligned} \frac{1}{r}\frac{\partial}{\partial \theta} \left( u_r {\boldsymbol{\mathbf{e}}}_r \right) \otimes {\boldsymbol{\mathbf{e}}}_{\theta} = \frac{1}{r}\frac{\partial u_r}{\partial \theta} {\boldsymbol{\mathbf{e}}}_r \otimes {\boldsymbol{\mathbf{e}}}_{\theta} + \frac{u_r}{r} \frac{\partial {\boldsymbol{\mathbf{e}}}_r}{\partial \theta} \otimes {\boldsymbol{\mathbf{e}}}_{\theta} = \frac{1}{r}\frac{\partial u_r}{\partial \theta} {\boldsymbol{\mathbf{e}}}_r \otimes {\boldsymbol{\mathbf{e}}}_{\theta} + \frac{u_r}{r} {\boldsymbol{\mathbf{e}}}_{\theta} \otimes {\boldsymbol{\mathbf{e}}}_{\theta}\end{aligned}\]
The final result, expressed as a tensor, is \[\label{eq:divu} \nabla{\boldsymbol{\mathbf{u}}}= \left[ \begin{array}{ccc} \frac{\partial u_r}{\partial r} & \frac{1}{r} \frac{\partial u_r}{\partial \theta} - \frac{u_{\theta}}{r} & \frac{\partial u_r}{\partial z} \\ \\ \frac{\partial u_{\theta}}{\partial r} & \frac{1}{r} \frac{\partial u_{\theta}}{\partial \theta} + \frac{u_{r}}{r} & \frac{\partial u_{\theta}}{\partial z} \\ \\ \frac{\partial u_z}{\partial r} & \frac{1}{r} \frac{\partial u_z}{\partial \theta} & \frac{\partial u_z}{\partial z} \end{array} \right]\]
Combining Eqs. (17), (18) and (20), we can finally write the components of the stress tensor \({\boldsymbol{\mathbf{\sigma}}}\): \[\begin{aligned} \sigma_{rr} =& -p + 2\mu \frac{\partial u_r}{\partial r},\\ \sigma_{\theta\theta} =& - p + 2\mu \left( \frac{1}{r} \frac{\partial u_{\theta}}{\partial \theta} + \frac{u_r}{r} \right),\\ \sigma_{zz} =& -p + 2\mu \frac{\partial u_z}{\partial z}, \\ \sigma_{r \theta} = \sigma_{\theta r} =& \mu \left( \frac{1}{r} \frac{\partial u_r}{\partial \theta} + \frac{\partial u_{\theta}}{\partial r} - \frac{u_{\theta}}{r}\right), \\ \sigma_{rz} = \sigma_{zr} = & \mu \left( \frac{\partial u_r}{\partial z} + \frac{\partial u_z}{\partial r}\right), \\ \sigma_{\theta z}=\sigma_{z \theta} = &\mu \left( \frac{1}{r} \frac{\partial u_z}{\partial \theta} + \frac{\partial u_{\theta}}{\partial z}\right).\end{aligned}\]
Now, substituting in the radial component of Eq. (14), we have:
\[\label{Eq:stress_strain_expanded} \begin{array}{rl} & \left( \frac{\sigma_{rr}}{r} + \frac{\partial \sigma_{rr}}{\partial r} + \frac{1}{r} \frac{\partial \sigma_{r\theta}}{\partial \theta} + \frac{\partial \sigma_{rz}}{\partial z} - \frac{\sigma_{\theta\theta}}{r} \right) \mathbf{e}_r \\ \\ &= \left[ \left( -\frac{p}{r} + \frac{2\mu}{r} \frac{\partial u_r}{\partial r} \right) + \left( -\frac{\partial p}{\partial r} + 2 \mu \frac{\partial^2 u_r}{\partial r^2} \right) \right. \\ \\ & +\left. \left( \frac{\mu}{r^2} \frac{\partial^2 u_r}{\partial \theta^2} + \frac{\mu}{r} \frac{\partial^2 u_{\theta}}{\partial r \partial \theta} - \frac{\mu }{r^2} \frac{\partial u_{\theta}}{\partial \theta} \right) \right. \\ \\ & + \left. \left( \mu \frac{\partial^2 u_r}{\partial z^2} + \mu \frac{\partial^2 u_z}{\partial r \partial z} \right) - \left( -\frac{p}{r} + \frac{2\mu}{r^2} \frac{\partial u_{\theta}}{\partial {\theta}} + 2\mu \frac{u_r}{r^2} \right) \right] \mathbf{e}_r \\ \\ &= -\frac{\partial p}{\partial r} + \mu \left[ \frac{2}{r} \frac{\partial u_r}{\partial r} + 2 \frac{\partial^2 u_r}{\partial r^2} + \frac{1}{r} \frac{\partial^2 u_{\theta}}{\partial r \partial \theta} \right. \\ \\ & \left. + \frac{1}{r^2} \frac{\partial^2 u_r}{\partial \theta^2} - \frac{1}{r^2} \frac{\partial u_{\theta}}{\partial \theta} + \left( \frac{\partial^2 u_r}{\partial z^2} + \frac{\partial^2 u_z}{\partial r \partial z} \right) - \frac{2}{r^2} \frac{\partial u_{\theta}}{\partial \theta} - \frac{2}{r^2} u_r \right] \mathbf{e}_r \end{array}\]
We rewrite now the first two differential terms in the square brackets in the following way: \[\begin{aligned} \frac{2}{r} \frac{\partial u_r}{\partial r} + 2 \frac{\partial^2 u_r}{\partial r^2} = \frac{1}{r}\frac{\partial}{\partial r} \left( r \frac{\partial u_r}{\partial r} \right) + \frac{1}{r} \frac{\partial u_r}{\partial r} + \frac{\partial^2 u_r}{\partial r^2} \\ = \frac{1}{r}\frac{\partial}{\partial r} \left( r \frac{\partial u_r}{\partial r} \right) + \frac{u_r}{r^2}+\frac{\partial}{\partial r} \left[ \frac{1}{r} \frac{\partial (r u_r)}{\partial r} \right] \end{aligned}\] and the third term as: \[\begin{aligned} \frac{1}{r}\frac{\partial^2 u_{\theta}}{\partial r \partial \theta} = \frac{\partial}{\partial r}\left( \frac{1}{r} \frac{\partial u_{\theta}}{\partial \theta}\right) + \frac{1}{r^2} \frac{\partial u_{\theta}}{\partial \theta}.\end{aligned}\] Now, substituting in Eq. (21) and collecting the terms, we can write the radial component as \[\begin{aligned} - \frac{\partial p}{\partial r} + \mu \left\{-\frac{u_r}{r^2}+ \frac{\partial}{\partial r} \left[ \frac{1}{r} \frac{\partial (r u_r)}{\partial r} + \frac{1}{r} \frac{\partial u_{\theta}}{\partial \theta}+ \frac{\partial u_z}{\partial z} \right] \right. \\ \left. + \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_r}{\partial r} \right) +\frac{1}{r^2}\frac{\partial^2 u_r}{\partial \theta^2} + \frac{\partial^2 u_r}{\partial z^2} - \frac{2}{r^2} \frac{\partial u_{\theta}}{\partial \theta} \right\}.\end{aligned}\]
For an incompressible fluid the term in the square brackets is equal to zero (continuity equation) and thus, combining the previous result with Eq. (9) we can write the following momentum equation for the radial component of the velocity: \[\begin{array}{c} \displaystyle \frac{\partial u_r}{\partial t} + u_r \frac{\partial u_r}{\partial r} + \frac{ u_{\theta}} {r} \frac{\partial u_r}{\partial \theta} - \frac{u_{\theta}^2}{r} + u_z \frac{\partial u_r}{\partial z} = \\ \\ \displaystyle - \frac{1}{\rho} \frac{\partial p}{\partial r} + \frac{\mu}{\rho} \left\{-\frac{u_r}{r^2} + \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_r}{\partial r} \right) +\frac{1}{r^2}\frac{\partial^2 u_r}{\partial \theta^2} + \frac{\partial^2 u_r}{\partial z^2} - \frac{2}{r^2} \frac{\partial u_{\theta}}{\partial \theta} \right\} + g_r. \end{array}\]
In a similar way, we can derive the momentum equations for the \(\theta\)-component of the velocity \[\begin{array}{c} \displaystyle \frac{\partial u_{\theta}}{\partial t} + u_r \frac{\partial u_{\theta}}{\partial r} + \frac{ u_{\theta}} {r} \frac{\partial u_{\theta}}{\partial \theta} + \frac{u_r u_{\theta}}{r} + u_z \frac{\partial u_{\theta}}{\partial z} = \\ \\ \displaystyle - \frac{1}{r\rho} \frac{\partial p}{\partial \theta} + \frac{\mu}{\rho} \left\{-\frac{u_{\theta}}{r^2} + \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_{\theta}}{\partial r} \right) +\frac{1}{r^2}\frac{\partial^2 u_{\theta}}{\partial \theta^2} + \frac{\partial^2 u_{\theta}}{\partial z^2} + \frac{2}{r^2} \frac{\partial u_r}{\partial \theta} \right\} + g_{\theta}, \end{array}\] and, for the \(z\)-component: \[\begin{array}{c} \displaystyle \frac{\partial u_z}{\partial t} + u_r \frac{\partial u_z}{\partial r} + \frac{ u_{\theta}} {r} \frac{\partial u_z}{\partial \theta} + u_z \frac{\partial u_z}{\partial z} = \\ \\ \displaystyle - \frac{1}{\rho} \frac{\partial p}{\partial z} + \frac{\mu}{\rho} \left\{ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_z}{\partial r} \right) +\frac{1}{r^2}\frac{\partial^2 u_z}{\partial \theta^2} + \frac{\partial^2 u_z}{\partial z^2} \right\} + g_z. \end{array}\]
Cylindrical coordinates are chosen to take advantage of symmetry, so that a velocity component can disappear. A very common case is axisymmetric flow with the assumption of no tangential velocity (\(u_{\theta}=0\)), and the remaining quantities are independent of \(\theta\). In this case the equations for the two remaining velocity components write as: \[\begin{array}{c} \displaystyle \frac{\partial u_r}{\partial t} + u_r \frac{\partial u_r}{\partial r} + u_z \frac{\partial u_r}{\partial z} = - \frac{1}{\rho} \frac{\partial p}{\partial r} + \frac{\mu}{\rho} \left\{-\frac{u_r}{r^2} + \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_r}{\partial r} \right) + \frac{\partial^2 u_r}{\partial z^2} \right\} + g_r, \end{array}\]
\[\begin{array}{c} \displaystyle \frac{\partial u_z}{\partial t} + u_r \frac{\partial u_z}{\partial r} + u_z \frac{\partial u_z}{\partial z} = - \frac{1}{\rho} \frac{\partial p}{\partial z} + \frac{\mu}{\rho} \left\{ \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_z}{\partial r} \right) + \frac{\partial^2 u_z}{\partial z^2} \right\} + g_z, \end{array}\]
while the continuity equation is \[\frac{1}{r}\frac{\partial}{\partial r} (r u_r) +\frac{\partial u_z}{\partial z} = 0.\]
The laminar flow through a pipe of uniform (circular) cross-section is known as Hagen–Poiseuille flow. The equations governing the Hagen–Poiseuille flow can be derived directly from the previous equations by making the following additional assumptions:
the flow is steady (\(\displaystyle \partial (...)/\partial t=0\));
the radial and swirl components of the velocity are zero (\(\displaystyle u_{r}=u_{\theta }=0\));
the flow is fully developed (\(\displaystyle \partial u_z / \partial z = 0\)).
The continuity equation is identically satisfied. The momentum equation for the radial component of the velocity reduces to \(\displaystyle \partial p/\partial r=0\), i.e., the pressure \(p\) is a function of the axial coordinate \(z\) only. The third momentum equation reduces to: \[\frac{1}{r}\frac{\partial}{\partial r} \left( r \frac{\partial u_z}{\partial r} \right) = \frac{1}{\mu} \frac{\partial p}{\partial z}.\] The equation can be integrated with respect to \(r\) and the solution is \[u_z = -\frac{1}{4\mu} \frac{\partial p}{\partial z} (R^2 - r^2 )\] where \(R\) is the radius of the pipe.
The maximum velocity occurs at the pipe centerline (\(r=0\)): \[\begin{aligned} u_{z,max} = \frac{R^2}{4\mu}\left( -\frac{\partial p}{\partial z} \right)\end{aligned}\] while the average velocity can be obtained by integrating over the pipe cross section and is \(u_{z,avg}= 0.5 u_{z,max}\).
From the velocity profile it is also possible to calculate the wall shear stress in the direction \(z\) (the wall shear stress in the direction \(\theta\) is zero, being \(u_r=u_{\theta}=0\)): \[\left. \sigma_{rz}\right|_{r=R}= \mu \left( \frac{\partial u_{r}}{\partial z} + \frac{\partial u_{z}}{\partial r} \right)_{r=R} = \mu \left( \frac{1}{4\mu} \frac{\partial p}{\partial z} 2r \right)_{r=R} = \frac{R}{2} \frac{\partial p}{\partial z}.\]
If we consider a pipe of length \(l\) and a constant vertical pressure gradient \(\frac{\partial p}{\partial z}=\Delta P/l\) (where \(\Delta P\) is the pressure difference between the outlet and inlet of the pipe), we can also calculate the shear stress at the wall from basic principles. The force resulting from the pressure difference is equal to \(\Delta p \pi R^2\) and should equal the opposing force generated by friction at the wall. This frictional force equals the shear stress \(\sigma_{rz}\), times the lateral surface \(2\pi R l\). Equating these forces gives \(\Delta P\cdot \pi R^2 = \sigma_{rz} \cdot 2\pi R l\), and thus \[\sigma_{rz} = \frac{R}{2} \frac{\partial p}{\partial z}.\]
In the divergence operator there is a factor \(1/r\) multiplying the partial derivative with respect to \(\theta\). An easy way to understand where this factor come from is to consider a function \(f(r,\theta,z)\) in cylindrical coordinates and its gradient. For a small change in going from a point \((r,\theta,z)\) to \((r+dr,\theta+d\theta,z+dz)\) we can write \[df = \frac{\partial f}{\partial r} dr + \frac{\partial f}{\partial \theta} d\theta + \frac{\partial f}{\partial z} dz.\] At the same time, we can express the same increment as \(df=(dr {\boldsymbol{\mathbf{e}}}_r + r d\theta {\boldsymbol{\mathbf{e}}}_{theta} + dz {\boldsymbol{\mathbf{e}}}_z)\cdot \nabla f\) Let assume that the gradient of f can be expressed as \[\nabla f = \alpha {\boldsymbol{\mathbf{e}}}_r + \beta {\boldsymbol{\mathbf{e}}}_{theta} + \gamma {\boldsymbol{\mathbf{e}}}_z\] where the scalar coefficients \(\alpha,\beta,\gamma\) are to be found. We have \[df = (dr {\boldsymbol{\mathbf{e}}}_r + r d\theta {\boldsymbol{\mathbf{e}}}_{theta} + dz {\boldsymbol{\mathbf{e}}}_z)\cdot ( \alpha {\boldsymbol{\mathbf{e}}}_r + \beta {\boldsymbol{\mathbf{e}}}_{theta} + \gamma {\boldsymbol{\mathbf{e}}}_z) = \alpha dr + r \beta d\theta + \gamma dz.\] Comparing this expression with the first one for \(df\) we see that \(\alpha=\frac{\partial f}{\partial r}\), \(\beta=\frac{1}{r}\frac{\partial f}{\partial \theta}\) and \(\gamma = \frac{\partial f}{\partial z}\) and thus we can write the gradient operator as \[\nabla = {\boldsymbol{\mathbf{e}}}_{r} \frac{\partial}{\partial r} + {\boldsymbol{\mathbf{e}}}_{\theta} \frac{1}{r}\frac{\partial}{\partial \theta} + {\boldsymbol{\mathbf{e}}}_{z} \frac{\partial}{\partial z}.\]↩