柱坐标系下的流体力学控制方程组的微分形式的推导

直角坐标系下描述

我们以NS方程（Navier-Stokes Equation）为例，来推导控制方程的柱坐标表示。我这里考虑不可压的，密度和粘性系数都为常数的情况。

此时，直角坐标系下的NS方程的表达形式为，
ρDV⃗Dt=−∇p+ρg⃗+μ∇2V⃗∇⋅V⃗=0\begin{array}{c}{\rho \frac{D \vec{V}}{D t}=-\nabla p+\rho \vec{g}+\mu \nabla^{2} \vec{V}} \\ {\nabla \cdot \vec{V}=0}\end{array} ρDtDV=−∇p+ρg+μ∇2V∇⋅V=0

按分量的形式可以写为：

ρ(∂u∂t+u∂u∂x+v∂u∂y+w∂u∂z)=−∂P∂x+ρgx+μ(∂2u∂x2+∂2u∂y2+∂2u∂z2)\rho\left(\frac{\partial u}{\partial t}+u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}+w \frac{\partial u}{\partial z}\right)=-\frac{\partial P}{\partial x}+\rho g_{x}+\mu\left(\frac{\partial^{2} u}{\partial x^{2}}+\frac{\partial^{2} u}{\partial y^{2}}+\frac{\partial^{2} u}{\partial z^{2}}\right) ρ(∂t∂u+u∂x∂u+v∂y∂u+w∂z∂u)=−∂x∂P+ρgx+μ(∂x2∂2u+∂y2∂2u+∂z2∂2u)
ρ(∂v∂t+u∂v∂x+v∂v∂y+w∂v∂z)=−∂P∂y+ρgy+μ(∂2v∂x2+∂2v∂y2+∂2v∂z2)\rho\left(\frac{\partial v}{\partial t}+u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+w \frac{\partial v}{\partial z}\right)=-\frac{\partial P}{\partial y}+\rho g_{y}+\mu\left(\frac{\partial^{2} v}{\partial x^{2}}+\frac{\partial^{2} v}{\partial y^{2}}+\frac{\partial^{2} v}{\partial z^{2}}\right) ρ(∂t∂v+u∂x∂v+v∂y∂v+w∂z∂v)=−∂y∂P+ρgy+μ(∂x2∂2v+∂y2∂2v+∂z2∂2v)
ρ(∂w∂t+u∂w∂x+v∂w∂y+w∂w∂z)=−∂P∂z+ρgz+μ(∂2w∂x2+∂2w∂y2+∂2w∂z2)\rho\left(\frac{\partial w}{\partial t}+u \frac{\partial w}{\partial x}+v \frac{\partial w}{\partial y}+w \frac{\partial w}{\partial z}\right)=-\frac{\partial P}{\partial z}+\rho g_{z}+\mu\left(\frac{\partial^{2} w}{\partial x^{2}}+\frac{\partial^{2} w}{\partial y^{2}}+\frac{\partial^{2} w}{\partial z^{2}}\right) ρ(∂t∂w+u∂x∂w+v∂y∂w+w∂z∂w)=−∂z∂P+ρgz+μ(∂x2∂2w+∂y2∂2w+∂z2∂2w)

直角坐标系到柱坐标系

如图所示，我们建立直角坐标系和柱坐标系。

由柱坐标系的定义，我们知道，直角坐标系和柱坐标系满足这样一个关系：

[rθz]=[x2+y2arctan⁡(y/x)z],0≤θ<2π\left[ \begin{array}{l}{r} \\ {\theta} \\ {z}\end{array}\right]=\left[ \begin{array}{c}{\sqrt{x^{2}+y^{2}}} \\ {\arctan (y / x)} \\ {z}\end{array}\right], \quad 0 \leq \theta<2 \pi ⎣⎡rθz⎦⎤=⎣⎡x2+y2arctan(y/x)z⎦⎤,0≤θ<2π

[xyz]=[rcos⁡θrsin⁡θz]\left[ \begin{array}{l}{x} \\ {y} \\ {z}\end{array}\right]=\left[ \begin{array}{c}{r \cos \theta} \\ {r \sin \theta} \\ {z}\end{array}\right] ⎣⎡xyz⎦⎤=⎣⎡rcosθrsinθz⎦⎤

假设在直角坐标系中，三个坐标轴方向的单位向量分别表示为i、j、k\mathbf{i、j、k}i、j、k，在柱坐标系中，三个轴向的单位向量表示为er,eθ,ez\mathbf{e_r,e_\theta,e_z}er,eθ,ez，那么，如图所示，我们将i、j、k\mathbf{i、j、k}i、j、k分解到柱坐标系中，将er,eθ,ez\mathbf{e_r,e_\theta,e_z}er,eθ,ez分解到直角坐标系中，可以得到二者之间的一个关系。

[ereθez]=[cos⁡θsin⁡θ0−sin⁡θcos⁡θ0001][ijk]\left[ \begin{array}{c}{\mathbf{e_r}} \\ {\mathbf{e_\theta}} \\ {\mathbf{e_z}}\end{array}\right]=\left[ \begin{array}{ccc}{\cos \theta} & {\operatorname{sin} \theta} & {0} \\ {-\sin \theta} & {\cos \theta} & {0} \\ {0} & {0} & {1}\end{array}\right] \left[ \begin{array}{l}{\mathbf{i}} \\ {\mathbf{j}} \\ {\mathbf{k}}\end{array}\right] ⎣⎡ereθez⎦⎤=⎣⎡cosθ−sinθ0sinθcosθ0001⎦⎤⎣⎡ijk⎦⎤

[ijk]=[cos⁡θ-sin⁡θ0sin⁡θcos⁡θ0001][ereθez]\left[ \begin{array}{l}{\mathbf{i}} \\ {\mathbf{j}} \\ {\mathbf{k}}\end{array}\right]=\left[ \begin{array}{ccc}{\cos \theta} & {\operatorname{-sin} \theta} & {0} \\ {\sin \theta} & {\cos \theta} & {0} \\ {0} & {0} & {1}\end{array}\right] \left[ \begin{array}{c}{\mathbf{e_r}} \\ {\mathbf{e_\theta}} \\ {\mathbf{e_z}}\end{array}\right] ⎣⎡ijk⎦⎤=⎣⎡cosθsinθ0-sinθcosθ0001⎦⎤⎣⎡ereθez⎦⎤

坐标系变换的雅克比矩阵体现了两个坐标系中变量的偏导关系，表示如下：

柱坐标系下常用算子表示

为了推导不可压NS方程在柱坐标下的表示形式，我们先推导一些常用算子的柱坐标表示（重要但后面不一定全会用到），即：
∇f=∂f∂rer+1r∂f∂θeθ+∂f∂zez∇⋅f=1r∂∂r(rfr)+1r∂fθ∂θ+∂fz∂z∇×f=(1r∂fz∂θ−∂fθ∂z)er+(∂fr∂z−∂fz∂r)eθ+1r(∂∂r(rfθ)−∂fr∂θ)ez∇2f=1r∂∂r(r∂f∂r)+1r2∂2f∂θ2+∂2f∂z2\begin{aligned} \nabla f &=\frac{\partial f}{\partial r} \mathbf{e_r}+\frac{1}{r} \frac{\partial f}{\partial \theta} \mathbf{e_\theta}+\frac{\partial f}{\partial z} \mathbf{e_z}\\ \nabla \cdot \boldsymbol{f} &=\frac{1}{r} \frac{\partial}{\partial r}\left(r f_{r}\right)+\frac{1}{r} \frac{\partial f_{\theta}}{\partial \theta}+\frac{\partial f_{z}}{\partial z} \\ \nabla \times \boldsymbol{f} &=\left(\frac{1}{r} \frac{\partial f_{z}}{\partial \theta}-\frac{\partial f_{\theta}}{\partial z}\right) \mathbf{e_r}+\left(\frac{\partial f_{r}}{\partial z}-\frac{\partial f_{z}}{\partial r}\right) \mathbf{e_\theta}+\frac{1}{r}\left(\frac{\partial}{\partial r}\left(r f_{\theta}\right)-\frac{\partial f_{r}}{\partial \theta}\right) \mathbf{e_z}\\ \nabla^{2} f &=\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial f}{\partial r}\right)+\frac{1}{r^{2}} \frac{\partial^{2} f}{\partial \theta^{2}}+\frac{\partial^{2} f}{\partial z^{2}} \end{aligned} ∇f∇⋅f∇×f∇2f=∂r∂fer+r1∂θ∂feθ+∂z∂fez=r1∂r∂(rfr)+r1∂θ∂fθ+∂z∂fz=(r1∂θ∂fz−∂z∂fθ)er+(∂z∂fr−∂r∂fz)eθ+r1(∂r∂(rfθ)−∂θ∂fr)ez=r1∂r∂(r∂r∂f)+r21∂θ2∂2f+∂z2∂2f

重要PS： 注意到，这里散度和旋度的定义用到的是的fff在柱坐标系下的分量来表示（右边），而一阶和二阶梯度中的fff是一个函数，这个函数既可以在直角坐标系下，也可以在柱坐标系下。如果f=V⃗\mathbf{f} = \vec Vf=V，每个分量表示速度沿三个方向的分量。

因为Tex公式敲起来挺费事的，我这里只证明一下第一个，其他类似。使用坐标系单位向量之间的关系和链式法则。
∇f=fx+fy+fz=[(frcosθ−fθ1rsinθ)cosθ+(frsinθ+fθ1rcosθ)sinθ]er+[(fr1rcosθ−fθ1rsinθ)−sinθ+(frsinθ+fθ1rcosθ)cosθ]eθ+fzez=frer+1rfθeθ+fzez\nabla f = f_x+f_y+f_z=[(f_r\mathrm{cos \theta-f_\theta\frac{1}{r}\mathrm{sin\theta}})\mathrm{cos \theta}+(f_r\mathrm{sin\theta+f_\theta\frac{1}{r}\mathrm{cos\theta}})\mathrm{sin\theta}]\mathbf{e_r}\\+[(f_r\frac{1}{r}\mathrm{cos \theta-f_\theta\frac{1}{r}\mathrm{sin\theta}})\mathrm{-sin\theta}+(f_r\mathrm{sin\theta+f_\theta\frac{1}{r}\mathrm{cos\theta}})\mathrm{cos\theta}]\mathbf{e_\theta}+f_z\mathbf{e_z}\\=f_r\mathbf{e_r}+\frac{1}{r}f_\theta \mathbf{e_\theta}+f_z\mathbf{e_z}∇f=fx+fy+fz=[(frcosθ−fθr1sinθ)cosθ+(frsinθ+fθr1cosθ)sinθ]er+[(frr1cosθ−fθr1sinθ)−sinθ+(frsinθ+fθr1cosθ)cosθ]eθ+fzez=frer+r1fθeθ+fzez

不可压NS方程柱坐标形式推导

有了上述的工具，我们下一步要做的将Navier-Stokes方程直角坐标表达中的各种算子替换成上述的极坐标表示，整理合并即可。
对于NS方程，对左端求物质导数，NS方程表达为：
ρ∂V⃗∂t+ρV⃗⋅∇V⃗=−∇p+ρg⃗+μ∇2V⃗{\rho \frac{\partial \vec{V}}{\partial t}+\rho\vec V \cdot \nabla \vec V=-\nabla p+\rho \vec{g}+\mu \nabla^{2} \vec{V}} ρ∂t∂V+ρV⋅∇V=−∇p+ρg+μ∇2V

外力项不涉及求导，不发生变换，只要替换相应的变量即可。所以我们只要计算时间导数项，对流项，压强项和粘性项，以及连续性方程。

主要思路

我们想做的事情是，在NS方程的分量表达式

当中，将u,v,wu,v,wu,v,w换成ur,uθ,uzu_r,u_\theta,u_zur,uθ,uz，并且我柱坐标系下的三个式子是关于er,eθ,ez\mathbf{e_r,e_\theta,e_z}er,eθ,ez三个方向的分量，即考虑将三个式子做线性组合，左边对左边，右边多右边：
[柱式1柱式2柱式3]=[cos⁡θsin⁡θ0−sin⁡θcos⁡θ0001][直式1直式2直式3]\left[ \begin{array}{c}{柱式1} \\ {柱式2} \\ {柱式3}\end{array}\right]=\left[ \begin{array}{ccc}{\cos \theta} & {\operatorname{sin} \theta} & {0} \\ {-\sin \theta} & {\cos \theta} & {0} \\ {0} & {0} & {1}\end{array}\right] \left[ \begin{array}{l}{直式1} \\ {直式2} \\ {直式3}\end{array}\right] ⎣⎡柱式1柱式2柱式3⎦⎤=⎣⎡cosθ−sinθ0sinθcosθ0001⎦⎤⎣⎡直式1直式2直式3⎦⎤

考虑直角坐标系向量到柱坐标系的替换：
V⃗=[uvw]=[cos⁡θ-sin⁡θ0sin⁡θcos⁡θ0001][uruθuz]\vec V = \left[ \begin{array}{l}{u} \\ {v} \\ {w}\end{array}\right]=\left[ \begin{array}{ccc}{\cos \theta} & {\operatorname{-sin} \theta} & {0} \\ {\sin \theta} & {\cos \theta} & {0} \\ {0} & {0} & {1}\end{array}\right] \left[ \begin{array}{c}{{u_r}} \\ {{u_\theta}} \\ {{u_z}}\end{array}\right] V=⎣⎡uvw⎦⎤=⎣⎡cosθsinθ0-sinθcosθ0001⎦⎤⎣⎡uruθuz⎦⎤

将其代入，并利用梯度和拉普拉斯算子在柱坐标下的表示公式（不全用到）：

∇f=∂f∂rer+1r∂f∂θeθ+∂f∂zez∇⋅f=1r∂∂r(rfr)+1r∂fθ∂θ+∂fz∂z∇×f=(1r∂fz∂θ−∂fθ∂z)er+(∂fr∂z−∂fz∂r)eθ+1r(∂∂r(rfθ)−∂fr∂θ)ez∇2f=1r∂∂r(r∂f∂r)+1r2∂2f∂θ2+∂2f∂z2\begin{aligned} \nabla f &=\frac{\partial f}{\partial r} \mathbf{e_r}+\frac{1}{r} \frac{\partial f}{\partial \theta} \mathbf{e_\theta}+\frac{\partial f}{\partial z} \mathbf{e_z}\\ \nabla \cdot \boldsymbol{f} &=\frac{1}{r} \frac{\partial}{\partial r}\left(r f_{r}\right)+\frac{1}{r} \frac{\partial f_{\theta}}{\partial \theta}+\frac{\partial f_{z}}{\partial z} \\ \nabla \times \boldsymbol{f} &=\left(\frac{1}{r} \frac{\partial f_{z}}{\partial \theta}-\frac{\partial f_{\theta}}{\partial z}\right) \mathbf{e_r}+\left(\frac{\partial f_{r}}{\partial z}-\frac{\partial f_{z}}{\partial r}\right) \mathbf{e_\theta}+\frac{1}{r}\left(\frac{\partial}{\partial r}\left(r f_{\theta}\right)-\frac{\partial f_{r}}{\partial \theta}\right) \mathbf{e_z}\\ \nabla^{2} f &=\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial f}{\partial r}\right)+\frac{1}{r^{2}} \frac{\partial^{2} f}{\partial \theta^{2}}+\frac{\partial^{2} f}{\partial z^{2}} \end{aligned} ∇f∇⋅f∇×f∇2f=∂r∂fer+r1∂θ∂feθ+∂z∂fez=r1∂r∂(rfr)+r1∂θ∂fθ+∂z∂fz=(r1∂θ∂fz−∂z∂fθ)er+(∂z∂fr−∂r∂fz)eθ+r1(∂r∂(rfθ)−∂θ∂fr)ez=r1∂r∂(r∂r∂f)+r21∂θ2∂2f+∂z2∂2f

另外链式法则需要用到两个坐标系变换的求导关系：

综上进行合并整理，即可。

连续性方程（质量守恒）

由柱坐标下的散度公式，可知不可压连续性条件为：
1r∂(rur)∂r+1r∂(uθ)∂θ+∂uz∂z=0\frac{1}{r} \frac{\partial\left(r u_{r}\right)}{\partial r}+\frac{1}{r} \frac{\partial\left(u_{\theta}\right)}{\partial \theta}+\frac{\partial u_{z}}{\partial z}=0 r1∂r∂(rur)+r1∂θ∂(uθ)+∂z∂uz=0

时间导数项

因为柱坐标变换是对空间的变换，不涉及时间变量，所以时间导数项在直角坐标系下不发生变化，做一个式分量的线性组合即可。
[∂ur∂t∂uθ∂t∂uz∂t]=[cos⁡θsin⁡θ0−sin⁡θcos⁡θ0001][∂u∂t∂v∂t∂w∂t]\left[ \begin{array}{c}{\frac{\partial u_r}{\partial t}} \\ {\frac{\partial u_\theta}{\partial t}} \\ {\frac{\partial u_z}{\partial t}}\end{array}\right]=\left[ \begin{array}{ccc}{\cos \theta} & {\operatorname{sin} \theta} & {0} \\ {-\sin \theta} & {\cos \theta} & {0} \\ {0} & {0} & {1}\end{array}\right] \left[ \begin{array}{l}{\frac{\partial u}{\partial t}} \\ {\frac{\partial v}{\partial t}} \\ {\frac{\partial w}{\partial t}}\end{array}\right] ⎣⎡∂t∂ur∂t∂uθ∂t∂uz⎦⎤=⎣⎡cosθ−sinθ0sinθcosθ0001⎦⎤⎣⎡∂t∂u∂t∂v∂t∂w⎦⎤

对流项

因为V⃗\vec VV是一个向量，我们可以分别考虑它的每一个分量：

ρV⃗⋅∇f=ρ(∂f∂ru+1r∂f∂θv+∂f∂zw)\rho\vec{V}\cdot \nabla f =\rho(\frac{\partial f}{\partial r} {u}+\frac{1}{r} \frac{\partial f}{\partial \theta} {v}+\frac{\partial f}{\partial z} {w})ρV⋅∇f=ρ(∂r∂fu+r1∂θ∂fv+∂z∂fw)

将fff分别替换为u,v,wu,v,wu,v,w，并将其替换为ur,uθ,uzu_r,u_\theta,u_zur,uθ,uz，合并化简，即可得到对流项。

ρ(ur∂ur∂r+uθr∂ur∂θ−uθ2r+uz∂ur∂z)\begin{array}{l}{\rho\left(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)} \end{array} ρ(ur∂r∂ur+ruθ∂θ∂ur−ruθ2+uz∂z∂ur)
ρ(ur∂uθ∂r+uθr∂uθ∂θ+uruθr+uz∂uθ∂z)\begin{array}{l}{\rho\left(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}\right)} \end{array} ρ(ur∂r∂uθ+ruθ∂θ∂uθ+ruruθ+uz∂z∂uθ)
ρ(ur∂uz∂r+uθr∂uz∂θ+uz∂uz∂z)\begin{array}{l}{\rho\left(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)} \end{array} ρ(ur∂r∂uz+ruθ∂θ∂uz+uz∂z∂uz)

打tex比较麻烦，将手算稿纸黏贴如下，(ρ\rhoρ在外面，先不考虑，下同)：

压强项

我们想要的其实就是直角坐标系中的各个分量在柱坐标系下的表达，由前提到的梯度公式，直接可得：

−∇p=−∂p∂rer−1r∂p∂θeθ−∂p∂zez-\nabla p =- \frac{\partial p}{\partial r} \mathbf{e_r}-\frac{1}{r} \frac{\partial p}{\partial \theta} \mathbf{e_\theta}-\frac{\partial p}{\partial z} \mathbf{e_z}−∇p=−∂r∂per−r1∂θ∂peθ−∂z∂pez

粘性项

∇2f=1r∂∂r(r∂f∂r)+1r2∂2f∂θ2+∂2f∂z2\nabla^{2} f =\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial f}{\partial r}\right)+\frac{1}{r^{2}} \frac{\partial^{2} f}{\partial \theta^{2}}+\frac{\partial^{2} f}{\partial z^{2}}∇2f=r1∂r∂(r∂r∂f)+r21∂θ2∂2f+∂z2∂2f

将fff分别替换为u,v,wu,v,wu,v,w，做式分量的线性组合，并使用ur,uθ,uzu_r,u_\theta,u_zur,uθ,uz线性表出替换后并整理，即可。
Tex打起来比较麻烦，手稿如下：

合并

因为等式两边貌似没有一个统一表示的方法，所以最后只能写成分量的形式，最后方程的形式可以写为：
ρ(∂ur∂t+ur∂ur∂r+uθr∂ur∂θ−uθ2r+uz∂ur∂z)=−∂P∂r+ρgr+μ[1r∂∂r(r∂ur∂r)−urr2+1r2∂2ur∂θ2−2r2∂uθ∂θ+∂2ur∂z2]\begin{array}{l}{\rho\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)} \\ {=-\frac{\partial P}{\partial r}+\rho g_{r}+\mu\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial u_{r}}{\partial r}\right)-\frac{u_{r}}{r^{2}}+\frac{1}{r^{2}} \frac{\partial^{2} u_{r}}{\partial \theta^{2}}-\frac{2}{r^{2}} \frac{\partial u_{\theta}}{\partial \theta}+\frac{\partial^{2} u_{r}}{\partial z^{2}}\right]}\end{array} ρ(∂t∂ur+ur∂r∂ur+ruθ∂θ∂ur−ruθ2+uz∂z∂ur)=−∂r∂P+ρgr+μ[r1∂r∂(r∂r∂ur)−r2ur+r21∂θ2∂2ur−r22∂θ∂uθ+∂z2∂2ur]

ρ(∂uθ∂t+ur∂uθ∂r+uθr∂uθ∂θ+uruθr+uz∂uθ∂z)=−1r∂P∂θ+ρgθ+μ[1r∂∂r(r∂uθ∂r)−uθr2+1r2∂2uθ∂θ2+2r2∂ur∂θ+∂2uθ∂z2]\begin{array}{l}{\rho\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_{r} u_{\theta}}{r}+u_{z} \frac{\partial u_{\theta}}{\partial z}\right)} \\ {=-\frac{1}{r} \frac{\partial P}{\partial \theta}+\rho g_{\theta}+\mu\left[\frac{1}{r} \frac{\partial}{\partial r}\left(r \frac{\partial u_{\theta}}{\partial r}\right)-\frac{u_{\theta}}{r^{2}}+\frac{1}{r^{2}} \frac{\partial^{2} u_{\theta}}{\partial \theta^{2}}+\frac{2}{r^{2}} \frac{\partial u_{r}}{\partial \theta}+\frac{\partial^{2} u_{\theta}}{\partial z^{2}}\right]}\end{array} ρ(∂t∂uθ+ur∂r∂uθ+ruθ∂θ∂uθ+ruruθ+uz∂z∂uθ)=−r1∂θ∂P+ρgθ+μ[r1∂r∂(r∂r∂uθ)−r2uθ+r21∂θ2∂2uθ+r22∂θ∂ur+∂z2∂2uθ]

ρ(∂uz∂t+ur∂uz∂r+uθr∂uz∂θ+uz∂uz∂z)=−∂P∂z+ρgz+μ[1r∂∂r(r∂uz∂r)+1r2∂2uz∂θ2+∂2uz∂z2]\begin{array}{l}{\rho\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)} \\ {=-\frac{\partial P}{\partial z}+\rho g_{z}+\mu\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]}\end{array} ρ(∂t∂uz+ur∂r∂uz+ruθ∂θ∂uz+uz∂z∂uz)=−∂z∂P+ρgz+μ[r1∂r∂(r∂r∂uz)+r21∂θ2∂2uz+∂z2∂2uz]

不可压连续性条件为：
1r∂(rur)∂r+1r∂(uθ)∂θ+∂uz∂z=0\frac{1}{r} \frac{\partial\left(r u_{r}\right)}{\partial r}+\frac{1}{r} \frac{\partial\left(u_{\theta}\right)}{\partial \theta}+\frac{\partial u_{z}}{\partial z}=0 r1∂r∂(rur)+r1∂θ∂(uθ)+∂z∂uz=0