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

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

直角坐标系下描述

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

此时，直角坐标系下的NS方程的表达形式为，
$\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​

按分量的形式可以写为：

$\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​)
$\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​)
$\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​)

直角坐标系到柱坐标系

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

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

$\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&lt;2 \pi$⎣⎡​rθz​⎦⎤​=⎣⎡​x2+y2 ​arctan(y/x)z​⎦⎤​,0≤θ<2π

$\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​⎦⎤​

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

$\left[ \begin{array}{c}{\mathbf{e_r}} \\ {\mathbf{e_\theta}} \\ {\mathbf{e_z}}\end{array}\right]=\left[ \begin{array}{ccc}{\cos \theta} &amp; {\operatorname{sin} \theta} &amp; {0} \\ {-\sin \theta} &amp; {\cos \theta} &amp; {0} \\ {0} &amp; {0} &amp; {1}\end{array}\right] \left[ \begin{array}{l}{\mathbf{i}} \\ {\mathbf{j}} \\ {\mathbf{k}}\end{array}\right]$⎣⎡​er​eθ​ez​​⎦⎤​=⎣⎡​cosθ−sinθ0​sinθcosθ0​001​⎦⎤​⎣⎡​ijk​⎦⎤​

$\left[ \begin{array}{l}{\mathbf{i}} \\ {\mathbf{j}} \\ {\mathbf{k}}\end{array}\right]=\left[ \begin{array}{ccc}{\cos \theta} &amp; {\operatorname{-sin} \theta} &amp; {0} \\ {\sin \theta} &amp; {\cos \theta} &amp; {0} \\ {0} &amp; {0} &amp; {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θ0​001​⎦⎤​⎣⎡​er​eθ​ez​​⎦⎤​

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

柱坐标系下常用算子表示

为了推导不可压NS方程在柱坐标下的表示形式，我们先推导一些常用算子的柱坐标表示（重要但后面不一定全会用到），即：
$\begin{aligned} \nabla f &amp;=\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} &amp;=\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} &amp;=\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 &amp;=\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∂f​er​+r1​∂θ∂f​eθ​+∂z∂f​ez​=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：**注意到，这里散度和旋度的定义用到的是的 $f$f在柱坐标系下的分量来表示，而一阶和二阶梯度中的 $f$f是一个函数，这个函数既可以在直角坐标系下，也可以在柱坐标系下。

因为Tex公式敲起来挺费事的，我这里只证明一下第一个，其他类似。使用坐标系单位向量之间的关系和链式法则。
$\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​=[(fr​cosθ−fθ​r1​sinθ)cosθ+(fr​sinθ+fθ​r1​cosθ)sinθ]er​+[(fr​r1​cosθ−fθ​r1​sinθ)−sinθ+(fr​sinθ+fθ​r1​cosθ)cosθ]eθ​+fz​ez​=fr​er​+r1​fθ​eθ​+fz​ez​

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

有了上述的工具，我们下一步要做的将Navier-Stokes方程直角坐标表达中的各种算子替换成上述的极坐标表示，整理合并即可。
对于NS方程，对左端求物质导数，NS方程表达为：
${\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方程的分量表达式

$\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​)
$\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​)
$\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​)

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

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

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

$\begin{aligned} \nabla f &amp;=\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} &amp;=\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} &amp;=\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 &amp;=\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∂f​er​+r1​∂θ∂f​eθ​+∂z∂f​ez​=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​​

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

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

连续性方程（质量守恒）

由柱坐标下的散度公式，可知不可压连续性条件为：
$\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

时间导数项

因为柱坐标变换是对空间的变换，不涉及时间变量，所以时间导数项在直角坐标系下不发生变化，做一个式分量的线性组合即可。
$\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} &amp; {\operatorname{sin} \theta} &amp; {0} \\ {-\sin \theta} &amp; {\cos \theta} &amp; {0} \\ {0} &amp; {0} &amp; {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θ0​sinθcosθ0​001​⎦⎤​⎣⎡​∂t∂u​∂t∂v​∂t∂w​​⎦⎤​

对流项

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

$\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∂f​u+r1​∂θ∂f​v+∂z∂f​w)

将 $f$f分别替换为 $u,v,w$u,v,w，并将其替换为 $u_r,u_\theta,u_z$ur​,uθ​,uz​，合并化简，即可得到对流项。

$\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​​)​
$\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θ​​+rur​​)​
$\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θ​​+rur​uθ​​+uz​∂z∂uθ​​)​)​
$\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$ρ在外面，先不考虑，下同)：

压强项

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

∂uz​​+