I'm trying to simulate a 2D axisymmetric model of steady-state compressible viscous flow using Mathematica, but I get some errors. There is a chance that I'm making some mistakes with the governing equations and/or the boundary conditions, as suggested in some posts (e.g. here).


The geometry is an axisymmetric tubular step about the x axis forming an annular orifice in 3D. enter image description here


  1. Axisymmetric about the x axis
  2. conduction and radiation are negligable
  3. Newtonian fluid with constant viscosity $\eta$
  4. Ideal gas (air)
  5. Steady-state

Writing the conservation equations in cylindrical form for the axisymmetric problem:

Conservation of mass:

$$ \frac{\partial}{\partial x}\left( \rho v_x \right)+\frac{1}{r}\frac{\partial }{\partial r}\left(r \rho v_r\right)=0 \tag{1}$$

Conservation of linear momentum in axial direction:

$$\frac{\partial}{\partial x}\left( \rho v_x^2+\mathring{R} \rho T \right)+\frac{1}{r}\frac{\partial}{\partial r}\left( r \left( \rho v_r v_x + \eta \frac{\partial v_x}{\partial r} \right)\right) \tag{2}$$

Conservation of linear momentum in radial direction:

$$ \frac{\partial}{\partial x}\left( \rho v_x v_r+\eta \frac{\partial v_r}{\partial x} \right)+ \frac{1}{r}\frac{\partial}{\partial r}\left( r \left( \rho v_r ^2 +\mathring{R} \rho T \right) \right)=0 \tag{3}$$

And energy balance (conservation of heat) more info here:

$$\rho c_v\left( v_x \frac{\partial T}{\partial x} + v_r \frac{\partial T}{\partial r} \right)+ \mathring{R} \rho T \left( \frac{1}{r}\frac{\partial}{\partial r} \left( r v_r \right)+ \frac{\partial v_x}{\partial x} \right)+ \eta \left( 2 \left( \frac{\partial v_x}{\partial x} \right)^2+ 2 \left( \frac{\partial v_r}{\partial r} \right)^2+ \left( \frac{\partial v_r}{\partial x}+ \frac{\partial v_x}{\partial r} \right)^2 \\ -\frac{2}{3}\left( \frac{1}{r} \frac{\partial}{\partial r}\left( r v_r \right) + \frac{\partial v_x}{\partial x} \right)^2 \right)=0 \tag{4}$$

The boundary conditions as I presume are:

  • (1) Inlet: uniform pressure $P=P_1$
  • (2) Outlet: uniform pressure $P=P_0$
  • (3) Axis of symmetry $v_r=0$ and $\frac{\partial *}{\partial r}=0$ for all variables $v_r$, $v_x$, $T$ and $\rho$
  • (4) no slip walls $v_*=0$

I would appreciate if you could help me know if I'm writing the equations correctly or I'm making any apparent mistakes? Are the boundary conditions enough/correct?

  • 1
    $\begingroup$ Minor comment on notation: $\nu$ is typically the kinematic viscosity ($\nu\equiv\mu/\rho$), so it's very confusing to see $\nu$ in place of $v$ in your equations. $\endgroup$
    – Kyle Kanos
    Feb 22, 2018 at 11:16
  • $\begingroup$ Well, Thanks. I edited the post. I will definitely use $v$ in future posts. $\endgroup$
    – Foad
    Feb 22, 2018 at 11:20
  • $\begingroup$ I'm not sure about some of your terms. is $\mathring{R}$ the universal gas constant? $\endgroup$
    – Derek
    Feb 22, 2018 at 16:24
  • $\begingroup$ Also, are you sure you need all the terms you are using? Do you really need the fully compressible Navier-Stokes, or is the Mach number low enough that treating the fluid as incompressible would be sufficient? Similarly, do you need the viscous heating terms? Have you estimated the Brinkman number to see if you can ignore that as well? $\endgroup$
    – Derek
    Feb 22, 2018 at 16:27
  • $\begingroup$ Can you walk me through your derivation a bit? I take it what you did is: 1. Take the continuity equation, N-S equations, and ideal gas law 2. Throw out all $\theta$ and $\frac{\partial}{\partial \theta}$ terms and equations. 3. Substitute all the $P$ terms in the continuity and N-S equations using the ideal gas law. Is that about right? $\endgroup$
    – Derek
    Feb 22, 2018 at 16:42

1 Answer 1


After doing quite some research I think I have the correct equations. Unfortunately all of the equations above (except for the continuity) are wrong.

Linear momentum:

The linear momentum conservation equation in the most compact and general form can be written as:

$$\rho \frac{D \check{v}}{D t}=\check{\nabla} \check{\sigma} \tag{5}$$

Where the $\check{\sigma}$ is the Cauchy stress tensor $\boldsymbol{\sigma}$ in matrix form. The Cauchy stress tensor can be divided into a volumetric/dilatational part and a deviatoric part:

$$\check{\sigma}= \check{\tau} -P\check{I} \tag{6}$$

In a perfect fluid P is the hydrostatic pressure and $\boldsymbol{\tau}$ is the viscous stress tensor.

Constitutive equation:

In an ideal Newtonian fluid viscous stress tensor can be calculated from:

$$ \check{\tau}=\eta\left( \check{\nabla}^T \check{v}+ \left( \check{\nabla}^T \check{v} \right)^T \right) +\lambda \left(\check{\nabla}\check{v}^T\right) \check{I} \tag{7}$$

Where $\eta$ is the absolute/dynamic viscosity and $\lambda$ is the bulk viscosity defined as:

$$\lambda=\kappa-\frac{2}{3}\eta \tag{8}$$

Where $\kappa$ is the dilation/expansion viscosity, which from Stokes’ hypothesis for monatomic gases at low density it is negligible impyling:

$$ \check{\tau}=\eta\left( \check{\nabla}^T \check{v}+ \left( \check{\nabla}^T \check{v} \right)^T -\frac{2}{3} \left(\check{\nabla}.\check{v}^T\right) \check{I} \right) \tag{9}$$

$\eta$ is not necessarily constant. There are some empirical models for example Sutherland:

$$\eta \approx C_S \frac{T^{\left(3/2\right)}}{T+T_S} \tag{10}$$

Where $C_S$ and $T_S$ are constants.

Combining the equations 5, 6 and 9 and also considering the axisymmetry, ideal gas and steady-state assumptions equations of linear momentum in axial and radial directions can be written as:

$$\begin{gather*} \rho \left( v_r \frac{\partial v_r}{\partial r} +v_x \frac{\partial v_r}{\partial x} \right) = \frac{\partial}{\partial r} \left( \eta \left( -\frac{2}{3}\left( \frac{1}{r}\frac{\partial}{\partial r} \left( r v_r \right) +\frac{\partial v_x}{\partial x} \right) +2 \frac{\partial v_r}{\partial r} \right) \right) \\ +\frac{\partial}{\partial x}\left( \eta\left( \frac{\partial v_r}{\partial x}+\frac{\partial v_x}{\partial r} \right) \right)+ \frac{2 \eta}{r}\left( \frac{\partial v_r}{\partial r}- \frac{v_r}{r} \right) -\frac{\partial }{\partial r}\left(\mathring{R}\rho T \right) \end{gather*} \tag{11}$$

In radial direction and

$$\begin{gather*} \rho \left( v_r \frac{\partial v_x}{\partial r} +v_x \frac{\partial v_x}{\partial x} \right) = \frac{\partial}{\partial x} \left( \eta \left( -\frac{2}{3}\left( \frac{1}{r}\frac{\partial}{\partial r} \left( r v_r \right) +\frac{\partial v_x}{\partial x} \right) +2 \frac{\partial v_x}{\partial x} \right) \right) \\ +\frac{\partial}{\partial r}\left( \eta\left( \frac{\partial v_r}{\partial x}+\frac{\partial v_x}{\partial r} \right) \right)+ \frac{ \eta}{r}\left( \frac{\partial v_r}{\partial x}+ \frac{\partial v_x}{\partial r} \right) -\frac{\partial }{\partial x}\left(\mathring{R}\rho T \right) \end{gather*} \tag{12}$$

In axial direction. (ref1, ref2)

Energy equation:

I'm still not completely sure about the correct form of the energy equation as I have explained here. But assuming my post here is valid then the energy equation in compact form can be written as:

$$ \rho \frac{D e}{D t}=\check{\sigma} : \check{\nabla}^T \check{v} \tag{13}$$

Considering that for an ideal gas $e=c_v T$ (more info here and here) and axisymmetry and steady-state assumptions it expands to:

$$ \begin{gather*} \rho c_v\left( v_x \frac{\partial T}{\partial x} + v_r \frac{\partial T}{\partial r} \right)+ \mathring{R} \rho T \left( \frac{1}{r}\frac{\partial}{\partial r} \left( r v_r \right)+ \frac{\partial v_x}{\partial x} \right)= \\ \eta \left( 2 \left( \frac{\partial v_x}{\partial x} \right)^2+ 2 \left( \frac{\partial v_r}{\partial r} \right)^2+ \left( \frac{\partial v_r}{\partial x}+ \frac{\partial v_x}{\partial r} \right)^2 -\frac{2}{3}\left( \frac{1}{r} \frac{\partial}{\partial r}\left( r v_r \right) + \frac{\partial v_x}{\partial x} \right)^2 \right) \end{gather*} \tag{14}$$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.