柱座標系下的流體力學控制方程組的微分形式的推導

柱座標系下的流體力學控制方程組的微分形式的推導

直角座標系下描述

我們以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​​+