柱坐标系下的N-S方程推导

标签：方程 frac 推导 cos sin rho theta 坐标系 partial

本文旨在采用最朴实的方式推导柱坐标系下的NS方程。 (x,y,z)⇒(r,θ,z)\left( x,y,z \right)\Rightarrow \left(r,\theta,z \right)

笛卡尔坐标系中，守恒形式的NS方程（不需要推导过程的话，直接拉到最下面看结果）：

连续性方程：

Eqn1 ∂ρ∂t+∇⋅(ρU)=0\frac{\partial \rho}{\partial t}+ \nabla \cdot \left( \rho \mathbf{U} \right) = 0

上式中时间项与坐标系无关，因此只需要对对流项进行转换：

Eqn2 ∇⋅(ρU)=∂ρUx∂x+∂ρUy∂y+∂ρUz∂z\nabla \cdot \left( \rho \mathbf{U} \right) = \frac{\partial \rho U_x}{\partial x} +\frac{\partial \rho U_y}{\partial y} +\frac{\partial \rho U_z}{\partial z}

先考虑等式右侧第一项。笛卡尔坐标系中x轴坐标与柱坐标的转换为 x=rcosθx = r cos \theta ，所以根据链式法则：

∂ρUx∂x=∂ρUx∂r∂r∂x+∂ρUx∂θ∂θ∂x\frac{\partial \rho U_x}{\partial x} = \frac{\partial \rho U_x}{\partial r}\frac{\partial r}{\partial x} + \frac{\partial \rho U_x}{\partial \theta} \frac{\partial \theta}{\partial x}

接下来从 x=rcosθx = r cos \theta 出发将上式中的四个偏导数推导出来：

1. ∂ρUx∂r \frac{\partial \rho U_x}{\partial r} 转换为柱坐标

DxDt=DrDtcosθ−rsinθDθDt⇒Ux=Urcosθ−sinθUθ\frac{Dx}{Dt} = \frac{Dr}{Dt} cos \theta - rsin\theta \frac{D\theta}{Dt}\Rightarrow U_x = U_r cos\theta - sin\theta U_\theta

⇒∂ρUx∂r=∂ρUr∂rcosθ+ρUr∂cosθ∂r−∂ρUθ∂rsinθ−ρUθ∂sinθ∂r\Rightarrow \frac{\partial \rho U_x}{\partial r} = \frac{\partial \rho U_r}{\partial r} cos \theta + \rho U_r \frac{\partial cos \theta}{\partial r} - \frac{\partial \rho U_\theta}{\partial r}sin\theta - \rho U_\theta\frac{\partial sin\theta}{\partial r}

又因为 ∂θ∂r=0\frac{\partial \theta}{\partial r}=0 ，因此上式简化为：

∂ρUx∂r=∂ρUr∂rcosθ−∂ρUθ∂rsinθ\frac{\partial \rho U_x}{\partial r} = \frac{\partial \rho U_r}{\partial r} cos \theta - \frac{\partial \rho U_\theta}{\partial r}sin\theta

2. ∂r∂x\frac{\partial r}{\partial x} 转换为柱坐标

r=x2+y2⇒∂r∂x=xx2+y2=rcosθ/r=cosθr=\sqrt{x^2+y^2} \Rightarrow \frac{\partial r}{\partial x} = \frac{x}{\sqrt{x^2+y^2}} =rcos\theta/r=cos\theta

3. ∂ρUx∂θ\frac{\partial \rho U_x}{\partial \theta} 转换为柱坐标

Ux=Urcosθ−sinθUθ⇒∂ρUx∂θ=∂ρUr∂θcosθ+ρUr∂cosθ∂θ−∂ρUθ∂θsinθ−ρUθsin⁡θ∂θU_x = U_r cos\theta - sin\theta U_\theta \Rightarrow \frac{\partial \rho U_x}{\partial \theta} = \frac{\partial \rho U_r}{\partial \theta} cos\theta + \rho U_r\frac{\partial cos\theta}{\partial \theta} - \frac{\partial \rho U_\theta}{\partial \theta} sin\theta - \rho U_\theta \frac{\sin\theta}{\partial \theta}

=∂ρUr∂θcosθ−ρUrsinθ−∂ρUθ∂θsinθ−ρUθcosθ=\frac{\partial \rho U_r}{\partial\theta}cos\theta - \rho U_r sin\theta - \frac{\partial \rho U_\theta}{\partial \theta} sin\theta - \rho U_\theta cos\theta

4. ∂θ∂x\frac{\partial \theta}{\partial x} 转换为柱坐标

θ=arctan(yx)⇒∂θ∂x=−yx2+y2=−rsinθr2=−sinθr\theta = arctan\left( \frac{y}{x} \right)\Rightarrow \frac{\partial \theta}{\partial x} = -\frac{y}{x^2+y^2} = -\frac{r sin\theta}{r^2} = -\frac{sin\theta}{r}

将1，2，3，4的结果全部带入式2

∂ρUx∂x=∂ρUx∂r∂r∂x+∂ρUx∂θ∂θ∂x\frac{\partial \rho U_x}{\partial x} = \frac{\partial \rho U_x}{\partial r}\frac{\partial r}{\partial x} + \frac{\partial \rho U_x}{\partial \theta} \frac{\partial \theta}{\partial x}

=∂ρUr∂rcos2θ−∂ρUθ∂rsinθcosθ−1r∂ρUr∂θcosθsinθ+1rρUrsin2θ+1r∂ρUθ∂θsin2θ+1rρUθcosθsinθ= \frac{\partial \rho U_r}{\partial r} cos^2 \theta - \frac{\partial \rho U_\theta}{\partial r}sin\theta cos\theta - \frac{1}{r} \frac{\partial \rho U_r}{\partial\theta}cos\theta sin\theta + \frac{1}{r} \rho U_r sin^2\theta + \frac{1}{r} \frac{\partial \rho U_\theta}{\partial \theta} sin^2\theta + \frac{1}{r} \rho U_\theta cos\theta sin\theta

同理可得：

∂ρUy∂y=∂ρUr∂rsin2θ+∂ρUθ∂rsinθcosθ+1r∂ρUr∂θsinθcosθ+1rρUrcos2θ+1r∂ρUθ∂θcos2θ−1rρUθsinθcosθ\frac{\partial \rho U_y}{\partial y} = \frac{\partial \rho U_r}{\partial r} sin^2\theta +\frac{\partial \rho U_\theta}{\partial r}sin\theta cos\theta + \frac{1}{r} \frac{\partial \rho U_r}{\partial \theta} sin\theta cos\theta + \frac{1}{r} \rho U_r cos^2\theta +\frac{1}{r} \frac{\partial \rho U_\theta}{\partial \theta} cos^2\theta -\frac{1}{r} \rho U_\theta sin\theta cos\theta

因此可得：

∂ρ∂t+∇⋅(ρU)=∂ρ∂t+∂ρUr∂r+1rρUr+1r∂ρUθ∂θ+∂ρUz∂z\frac{\partial \rho}{\partial t}+ \nabla \cdot \left( \rho \mathbf{U} \right) = \frac{\partial \rho}{\partial t} + \frac{\partial \rho U_r}{\partial r} +\frac{1}{r} \rho U_r +\frac{1}{r} \frac{\partial \rho U_\theta}{\partial \theta} +\frac{\partial \rho U_z}{\partial z}

动量方程

∂ρU∂t+∇⋅(ρUU)=−∇p+∇⋅τ+F\frac{\partial \rho \mathbf{U}}{\partial t} + \nabla \cdot \left( \rho \mathbf{U} \mathbf{U} \right) = -\nabla p +\nabla\cdot \mathbf{\tau}+\mathbf{F}

方程左侧：

其中p为压力， τ\tau 为应力张量。先看方程左侧。对流项中将三个方向拆开写：

x方向：

Eqn3 ∂ρux∂t+∂ρuxux∂x+∂ρuxuy∂y+∂ρuxuz∂z\frac{\partial \rho u_x}{\partial t} + \frac{\partial \rho u_x u_x}{\partial x} +\frac{\partial \rho u_x u_y}{\partial y} + \frac{\partial \rho u_x u_z}{\partial z}

y方向：

Eqn4 ∂ρuy∂t+∂ρuyux∂x+∂ρuyuy∂y+∂ρuyuz∂z\frac{\partial \rho u_y}{\partial t} +\frac{\partial \rho u_y u_x}{\partial x} +\frac{\partial \rho u_y u_y}{\partial y} +\frac{\partial \rho u_y u_z}{\partial z}

由于 ur=cosθux+sinθuyu_r = cos\theta u_x + sin\theta u_y 因此将 cosθ⋅Eqn3+sinθ⋅Eqn4cos\theta\cdot Eqn3+sin\theta\cdot Eqn4 可得：

∂ρuxcosθ∂t=−ρuxsinθ∂θ∂t+cosθ∂ρux∂t\frac{\partial \rho u_x cos\theta}{\partial t} = -\rho u_xsin\theta\frac{\partial \theta}{\partial t} +cos\theta \frac{\partial \rho u_x}{\partial t}

⇒\Rightarrow cosθ∂ρux∂t=∂ρuxcosθ∂t+ρuxsinθ∂θ∂tcos\theta \frac{\partial \rho u_x}{\partial t} =\frac{\partial \rho u_x cos\theta}{\partial t} + \rho u_x sin\theta \frac{\partial \theta}{\partial t}

cosθ⋅Eqn3+sinθ⋅Eqn4=∂ρ(uxcosθ+uysinθ)∂t+∂ρ(uxcosθ+uysinθ)ux∂x+∂ρ(uxcosθ+uysinθ)uy∂y+∂ρ(uxcosθ+uysinθ)uz∂z+ρ(uxsinθ−uycosθ)(∂θ∂t+ux∂θ∂x+uy∂θ∂y+uz∂θ∂z)cos\theta\cdot Eqn3+sin\theta\cdot Eqn4\\ = \frac{\partial \rho (u_x cos\theta + u_y sin\theta)}{\partial t} +\frac{\partial \rho (u_x cos\theta + u_y sin\theta) u_x}{\partial x} + \frac{\partial \rho (u_x cos\theta + u_y sin\theta) u_y}{\partial y} +\frac{\partial \rho (u_x cos\theta + u_y sin\theta) u_z}{\partial z} \\ + \rho (u_x sin\theta - u_y cos\theta) \left( \frac{\partial \theta}{\partial t} +u_x\frac{\partial \theta}{\partial x} +u_y\frac{\partial \theta}{\partial y} +u_z\frac{\partial\theta}{\partial z} \right)

=∂ρur∂t+∇⋅(ρuru)−ρuθDθDt=\frac{\partial \rho u_r}{\partial t} + \nabla \cdot \left( \rho u_r \mathbf{u} \right) -\rho u_\theta\frac{D\theta}{Dt}

上式中前两项与连续性方程一样，只不过是标量从连续性方程中的 ρ\rho 变成了 ρur\rho u_r

所以最终有：

r方向： ∂ρur∂t+∂ρurur∂r+ρururr+∂ρuruθr∂θ+∂ρuruz∂z−ρuθ2r\frac{\partial \rho u_r}{\partial t} + \frac{\partial \rho u_r u_r}{\partial r} + \frac{\rho u_r u_r}{r} + \frac{\partial \rho u_r u_\theta}{r\partial \theta} + \frac{\partial \rho u_r u_z}{\partial z} - \frac{\rho u_\theta^2}{r}

同理可得 θ\theta 方向上用 −sinθ⋅Eqn3+cosθ⋅Eqn4-sin\theta \cdot Eqn3 + cos\theta \cdot Eqn4 去构建：

∂ρuθ∂t+∇⋅(ρuθu)−ρuruθr=∂ρuθ∂t+∂ρuruθ∂r+1rρuruθ+∂ρuθuθ∂θ+∂uθuz∂z+ρuruθr\frac{\partial \rho u_\theta}{\partial t} + \nabla\cdot (\rho u_\theta \mathbf{u}) - \frac{\rho u_r u_\theta}{r} \\ =\frac{\partial \rho u_\theta}{\partial t} +\frac{\partial \rho u_r u_\theta}{\partial r} + \frac{1}{r}\rho u_r u_\theta +\frac{\partial \rho u_\theta u_\theta}{\partial \theta} +\frac{\partial u_\theta u_z}{\partial z} +\frac{\rho u_r u_\theta}{r}

z方向没有变化就不写了。

压力梯度项

用相同的思路：

r方向： cosθ∂p∂x+sinθ∂p∂y=cosθ(∂p∂r∂r∂x+∂p∂θ∂θ∂x)+sinθ(∂p∂r∂r∂y+∂p∂θ∂θ∂y)=∂p∂rcos\theta \frac{\partial p}{\partial x} + sin\theta\frac{\partial p}{\partial y} = cos \theta \left( \frac{\partial p}{\partial r} \frac{\partial r}{\partial x} + \frac{\partial p}{\partial \theta} \frac{\partial \theta}{\partial x} \right) + sin\theta \left( \frac{\partial p}{\partial r} \frac{\partial r}{\partial y} +\frac{\partial p}{\partial \theta} \frac{\partial \theta}{\partial y} \right) \\ =\frac{\partial p}{\partial r}

θ\theta 方向： ∂pr∂θ\frac{\partial p}{r\partial \theta}

应力张量太繁琐了，就懒得写了。最终结果如下：

结果：

连续性方程：

∂ρ∂t+∇⋅(ρU)=∂ρ∂t+∂ρUr∂r+1rρUr+1r∂ρUθ∂θ+∂ρUz∂z\frac{\partial \rho}{\partial t}+ \nabla \cdot \left( \rho \mathbf{U} \right) = \frac{\partial \rho}{\partial t} + \frac{\partial \rho U_r}{\partial r} +\frac{1}{r} \rho U_r +\frac{1}{r} \frac{\partial \rho U_\theta}{\partial \theta} +\frac{\partial \rho U_z}{\partial z}

动量方程：

r方向： ∂ρur∂t+∂ρurur∂r+ρururr+∂ρuruθr∂θ+∂ρuruz∂z−ρuθ2r=−∂p∂r+μ[∂∂r(1r∂rur∂r+1r2∂2ur∂θ2+∂2ur∂z2−2r2∂uθ∂θ)]+fr\frac{\partial \rho u_r}{\partial t} + \frac{\partial \rho u_r u_r}{\partial r} + \frac{\rho u_r u_r}{r} + \frac{\partial \rho u_r u_\theta}{r\partial \theta} + \frac{\partial \rho u_r u_z}{\partial z} - \frac{\rho u_\theta^2}{r} \\=-\frac{\partial p}{\partial r} +\mu \left[ \frac{\partial }{\partial r} \left( \frac{1}{r} \frac{\partial r u_r}{\partial r} +\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) \right] + f_r

θ\theta 方向：

∂ρuθ∂t+∂ρuruθ∂r+1rρuruθ+∂ρuθuθ∂θ+∂uθuz∂z+ρuruθr=−1r∂p∂θ+μ[∂∂r(∂ruθr∂r)+∂2uθr2∂θ2+∂2uθ∂z2+2∂urr2∂θ]+fθ\frac{\partial \rho u_\theta}{\partial t} +\frac{\partial \rho u_r u_\theta}{\partial r} + \frac{1}{r}\rho u_r u_\theta +\frac{\partial \rho u_\theta u_\theta}{\partial \theta} +\frac{\partial u_\theta u_z}{\partial z} +\frac{\rho u_r u_\theta}{r} \\ =-\frac{1}{r}\frac{\partial p}{\partial \theta} +\mu \left[ \frac{\partial }{\partial r} \left( \frac{\partial r u_\theta}{r\partial r} \right) +\frac{\partial^2 u_\theta}{r^2\partial \theta^2} +\frac{\partial ^2 u_\theta}{\partial z^2} +\frac{2\partial u_r}{r^2\partial \theta} \right] + f_\theta

z方向：

∂ρuz∂t+∂ρuruz∂r+∂ρuθuzr∂θ+∂ρuzuz∂z=−∂p∂z+μ[1r∂∂r(r∂uz∂r)+∂2uzr2∂θ2+∂2uz∂z2]+fz\frac{\partial \rho u_z}{\partial t} +\frac{\partial \rho u_r u_z}{\partial r} +\frac{\partial \rho u_\theta u_z}{r\partial\theta} +\frac{\partial \rho u_z u_z}{\partial z} = -\frac{\partial p}{\partial z} + \mu\left[ \frac{1}{r}\frac{\partial}{\partial r} \left( r\frac{\partial u_z}{\partial r} \right) +\frac{\partial^2 u_z}{r^2\partial \theta^2} +\frac{\partial^2 u_z}{\partial z^2} \right] + f_z

 

标签：方程,frac,推导,cos,sin,rho,theta,坐标系,partial
From： https://www.cnblogs.com/anzhili/p/16879589.html