円柱座標系での連続体の運動方程式

先月、円柱座標系での連続体の運動方程式の導出を扱ったのですが、複数ある導出方法の中で初学者にとって理解しやすい方法についての解説がどこにも載っていなかったので自分で書くことにしました。

はじめに

本記事ではベクトルと行列のそれぞれを明示するために行列Aの場合は$[A]$、ベクトルAの場合は$\{A\}$のように[]{}を使い分けています。

また、本記事では円柱座標系での連続体の運動方程式の導出について簡単な紹介しか行わないため、詳細については下記の文献などを参照してください。

 弾性力学

amazon.co.jp

なお、この記事を読むにあたって以前書いたオイラーの運動方程式関連の話応力テンソルの内容を理解しておいた方が理解しやすいかもしれません。

また、ここでいう円柱座標系は下の図のような$r, \theta, z$座標系を指しています。 円柱座標系のポンチ絵

直交座標系での運動方程式

円筒座標系での運動方程式の導出方法は幾つかありますが、本記事では直交座標系での運動方程式を座標変換する方法を用います。そこでまず、次の図のような微小体積に働く応力と慣性力のつり合いを考えます。なお、簡単のため体積力は無視します。

微小体積

密度$\rho$は一定であるとし、$x$軸方向の変位を$u$とすると、$x$軸方向の運動方程式は

\begin{eqnarray} \rho dx dy dz\frac{\partial^2 u}{\partial t^2} &=& -\sigma_x dy dz + (\sigma_x + \frac{\partial \sigma_x}{\partial x} dx)dy dz - \tau_{xy} dx dz + (\tau_{xy} + \frac{\partial \tau_{xy}}{\partial y} dy)dx dz- \tau_{zx} dy dx + (\tau_{zx} + \frac{\partial \tau_{zx}}{\partial z}dz)dy dx \nonumber\\ \rho\frac{\partial^2 u}{\partial t^2} &=& \frac{\partial \sigma_x}{\partial x} + \frac{\partial \tau_{xy}}{\partial y} + \frac{\partial \tau_{zx}}{\partial z} \end{eqnarray}

となります。さらに、$y, z$軸方向の変位を$v, w$とし、先ほどと同様にそれぞれの方向の運動方程式をまとめると次の通りとなります。

\begin{eqnarray} \rho\frac{\partial^2 v}{\partial t^2} &=& \frac{\partial \sigma_y}{\partial y} + \frac{\partial \tau_{xy}}{\partial x} + \frac{\partial \tau_{yz}}{\partial z} \\ \rho\frac{\partial^2 w}{\partial t^2} &=& \frac{\partial \sigma_z}{\partial z} + \frac{\partial \tau_{zx}}{\partial x} + \frac{\partial \tau_{yz}}{\partial y} \end{eqnarray}

当たり前ですが、コーシーの運動方程式から体積力を無視したものとなります。

円柱座標系での運動方程式の導出

応力テンソルと慣性力の座標変換

詳細は文献[1]1に載っているので詳細は省略しますが、直交座標系から円柱座標系への座標変換行列$[P]$は次の通り表せます。

\begin{equation} [P] = \begin{bmatrix} \cos{\theta}&\sin{\theta}&0 \\ -\sin{\theta}&\cos{\theta}&0 \\ 0 & 0 &1 \end{bmatrix} \end{equation}

なお、この行列$[P]$は単純に$z$軸周りに角度$+\theta $だけ回転させた際の回転行列と同じです。$[P]$を用いて、座標変換前と変換後の応力テンソル$[\sigma]$と$[\sigma’]$の関係を表します。

\begin{equation} [\sigma] = \begin{bmatrix} \sigma_{x} & \tau_{xy} & \tau_{xz} \\ \tau_{yx} & \sigma_{y} & \tau_{yz} \\ \tau_{zx} & \tau_{zy} & \sigma_{z} \end{bmatrix} \end{equation}

\begin{equation} [\sigma^\prime]= \begin{bmatrix} \sigma_{r}&\tau_{r\theta}&\tau_{rz}\\ \tau_{\theta r}&\sigma_{\theta}&\tau_{\theta z}\\ \tau_{zr}&\tau_{z\theta}&\sigma_{z} \end{bmatrix} \end{equation}

任意の二階のテンソルについて、座標変換前のテンソル$[T]$と変換後のテンソル$[T’]$の関係は座標変換行列$[P]$を用いて次の通り表せます2

\begin{equation} [T’]=[P][T][P]^T \end{equation}

当然ながら式(7)は応力テンソルにおいても成り立ち、座標変換行列$[P]$は$[P]^{-1}=[P]^T$が成り立つため、

\begin{eqnarray} [\sigma’]&=&[P][\sigma][P]^T\nonumber\\ [\sigma]&=&[P]^T[\sigma’][P]\nonumber\\ &=& \begin{bmatrix} \cos{\theta}&-\sin{\theta}&0\\ \sin{\theta}&\cos{\theta}&0\\ 0&0&1 \end{bmatrix} \begin{bmatrix} \sigma_{r}&\tau_{r\theta}&\tau_{rz}\\ \tau_{\theta r}&\sigma_{\theta}&\tau_{\theta z}\\ \tau_{zr}&\tau_{z\theta}&\sigma_{z} \end{bmatrix} \begin{bmatrix} \cos{\theta}&\sin{\theta}&0\\ -\sin{\theta}&\cos{\theta}&0\\ 0&0&1 \end{bmatrix} \nonumber\\ &=& \begin{bmatrix} \sigma_{r}\cos^2{\theta}-\tau_{\theta r}\sin{2\theta}+\sigma_{\theta}\sin^2{\theta}&(\sigma_{r}-\sigma_{\theta})\sin{\theta}\cos{\theta}+\tau_{\theta r}\cos{2\theta}&\tau_{rz}\cos{\theta}-\tau_{\theta z}\sin{\theta}\\ (\sigma_{r}-\sigma_{\theta})\sin{\theta}\cos{\theta}+\tau_{\theta r}\cos{2\theta}&\sigma_{r}\sin^2{\theta}+\tau_{\theta r}\sin{2\theta}+\sigma_{\theta}\cos^2{\theta}&\tau_{rz}\sin{\theta}+\tau_{\theta z}\cos{\theta}\\ \tau_{rz}\cos{\theta}-\tau_{\theta z}\sin{\theta}&\tau_{rz}\sin{\theta}+\tau_{\theta z}\cos{\theta}&\sigma_{z} \end{bmatrix} \end{eqnarray}

となります。なお、式(8)を導出する際に応力テンソルの対称性は保たれているとしており、次の式を断ることなく用いています。

\begin{equation} \tau_{r\theta}=\tau_{\theta r},\ \ \tau_{rz}=\tau_{zr},\ \ \tau_{r\theta}=\tau_{\theta r},\ \ \tau_{rz}=\tau_{zr} \end{equation}

式(8)では行列の成分の個数としては9個存在しますが、応力テンソルの対称性から独立な成分は6個となります。したがって式(8)は

\begin{equation} \begin{Bmatrix} \sigma_x \\ \sigma_y \\ \sigma_z \\ \tau_{xy} \\ \tau_{xz} \\ \tau_{yz} \end{Bmatrix} = \begin{Bmatrix} \sigma_{r}\cos^2{\theta}-\tau_{\theta r}\sin{2\theta}+\sigma_{\theta}\sin^2{\theta} \\ \sigma_{r}\sin^2{\theta}+\tau_{\theta r}\sin{2\theta}+\sigma_{\theta}\cos^2{\theta} \\ \sigma_{z} \\ \left(\sigma_{r}-\sigma_{\theta}\right)\sin{\theta}\cos{\theta}+\tau_{\theta r}\cos{2\theta} \\ \tau_{rz}\cos{\theta}-\tau_{\theta z}\sin{\theta} \\ \tau_{rz}\sin{\theta}+\tau_{\theta z}\cos{\theta} \end{Bmatrix} \end{equation}

となります。続いて慣性力も座標変換を行います。円柱座標系での慣性力のベクトルを$\{F_{r,\ \theta,\ z}\}$とし、次の通り定めます。

\begin{equation} \begin{Bmatrix} F_{r,\ \theta,\ z} \end{Bmatrix}= \begin{Bmatrix} \rho\frac{\partial^2u_r}{\partial t^2} \\ \rho\frac{\partial^2u_\theta}{\partial t^2} \\ \rho\frac{\partial^2w}{\partial t^2} \end{Bmatrix} \end{equation}

直交座標系での慣性力を$\{F_{x,\ y,\ \ z}\}$としてベクトルとして表すと、

\begin{equation} \begin{Bmatrix} F_{x,\ y,\ \ z} \end{Bmatrix}= \begin{Bmatrix} \rho\frac{\partial^2u}{\partial t^2} \\ \rho\frac{\partial^2v}{\partial t^2} \\ \rho\frac{\partial^2w}{\partial t^2} \end{Bmatrix} \end{equation}

となります。これらの慣性力は単なるベクトルなので、両者の関係は

\begin{eqnarray} \{F_{x,\ y,\ \ z}\}&=& [P]^T\{F_{r,\ \theta,\ z}\} \nonumber\\ &=& \begin{bmatrix} \cos{\theta}&-\sin{\theta}&0\\ \sin{\theta}&\cos{\theta}&0\\ 0&0&1 \end{bmatrix} \begin{Bmatrix} \rho\frac{\partial^2u_r}{\partial t^2} \\ \rho\frac{\partial^2u_\theta}{\partial t^2} \\ \rho\frac{\partial^2w}{\partial t^2} \end{Bmatrix}= \begin{Bmatrix} \rho\frac{\partial^2u_r}{\partial t^2}\cos{\theta}-\rho\frac{\partial^2u_\theta}{\partial t^2}\sin{\theta} \\ \rho\frac{\partial^2u_r}{\partial t^2}\sin{\theta}+\rho\frac{\partial^2u_\theta}{\partial t^2}\cos{\theta} \\ \rho\frac{\partial^2w}{\partial t^2} \end{Bmatrix} \end{eqnarray}

で表せます。

偏微分演算子の変換

直交座標系での偏微分演算子を円柱座標上で表します。まず、下記の関係を導きます。

\begin{equation} \left\{ \begin{matrix} \frac{\partial r}{\partial x}=\frac{\partial}{\partial x}\left(\sqrt{x^2+y^2+z^2}\right)=\frac{2x}{2\sqrt{x^2+y^2+z^2}}=\frac{x}{r}=\frac{r\cos{\theta}}{r}=\cos{\theta}\\ \frac{\partial r}{\partial y}=\frac{\partial}{\partial y}\left(\sqrt{x^2+y^2+z^2}\right)=\frac{2y}{2\sqrt{x^2+y^2+z^2}}=\frac{y}{r}=\frac{r\sin{\theta}}{r}=\sin{\theta}\\ \frac{\partial r}{\partial z}=\frac{\partial}{\partial y}\left(\sqrt{x^2+y^2+z^2}\right)=\frac{2z}{2\sqrt{x^2+y^2+z^2}}=\frac{z}{r}=0\\ \frac{\partial\theta}{\partial x}=\frac{\partial\left(\tan^{-1}{\frac{y}{x}}\right)}{\partial x}=\frac{1}{\left(\frac{y}{x}\right)^2+1}\cdot\frac{d\left(\frac{y}{x}\right)}{dx}=\frac{x^2}{y^2+x^2}\cdot y\cdot\left(-\frac{1}{x^2}\right)=\frac{-r\sin{\theta}}{r^2}=\frac{-\sin{\theta}}{r}\\ \frac{\partial\theta}{\partial y}=\frac{\partial\left(\tan^{-1}{\frac{y}{x}}\right)}{\partial y}=\frac{1}{\left(\frac{y}{x}\right)^2+1}\cdot\frac{d\left(\frac{y}{x}\right)}{dy}=\frac{x^2}{y^2+x^2}\cdot\frac{1}{x}=\frac{r\cos{\theta}}{r^2}=\frac{\cos{\theta}}{r}\\ \frac{\partial\theta}{\partial z}=\frac{\partial z}{\partial\theta}=0\\ \frac{\partial z}{\partial z}=\frac{\partial z}{\partial z} \end{matrix} \right. \end{equation}

式(14)を用いて、直交座標系での偏微分演算子$(\frac{\partial}{\partial x},\ \frac{\partial}{\partial y},\ \frac{\partial}{\partial z})$を変形します。

\begin{equation} \left\{ \begin{matrix} \frac{\partial}{\partial x}=\frac{\partial r}{\partial x}\frac{\partial}{\partial r}+\frac{\partial\theta}{\partial x}\frac{\partial}{\partial\theta}+\frac{\partial z}{\partial x}\frac{\partial}{\partial z}=\cos{\theta}\frac{\partial}{\partial r}-\frac{\sin{\theta}}{r}\frac{\partial}{\partial\theta} \\ \frac{\partial}{\partial y}=\frac{\partial r}{\partial y}\frac{\partial}{\partial r}+\frac{\partial\theta}{\partial y}\frac{\partial}{\partial\theta}+\frac{\partial z}{\partial y}\frac{\partial}{\partial z}=\sin{\theta}\frac{\partial}{\partial r}+\frac{\cos{\theta}}{r}\frac{\partial}{\partial\theta}\\ \frac{\partial}{\partial z}=\frac{\partial}{\partial z} \end{matrix} \right. \end{equation}

運動方程式の導出

直交座標系での$x$軸方向の運動方程式である式(1)に式(10), (13), (15)を代入して整理します。

\begin{eqnarray} \rho\frac{\partial^2 u_r}{\partial t^2}\cos{\theta}-\rho\frac{\partial^2u_\theta}{\partial t^2}\sin{\theta}&=& \left(\cos{\theta}\frac{\partial}{\partial r}-\frac{\sin{\theta}}{r}\frac{\partial}{\partial\theta}\right) \left(\sigma_{r}\cos^2{\theta}-\tau_{\theta r}\sin{2\theta}+\sigma_{\theta}\sin^2{\theta}\right)+ \left(\sin{\theta}\frac{\partial}{\partial r}+\frac{\cos{\theta}}{r}\frac{\partial}{\partial\theta}\right)\left\{\left(\sigma_{r}-\sigma_{\theta}\right)\sin{\theta}\cos{\theta}+\tau_{\theta r}\cos{2\theta}\right\}+\frac{\partial}{\partial z} \left(\tau_{rz}\cos{\theta}-\tau_{\theta z}\sin{\theta}\right) \nonumber\\ \frac{\partial^2u_r}{\partial t^2}\cos{\theta}-\rho\frac{\partial^2u_\theta}{\partial t^2}\sin{\theta}&=& \frac{\partial\sigma_{r}}{\partial r}\cos{\theta}-\frac{\partial\tau_{\theta r}}{\partial r}\sin{\theta}+\frac{1}{r}\left(\frac{\partial\tau_{\theta r}}{\partial\theta}\cos{\theta}-\frac{\partial\sigma_{\theta}}{\partial\theta}\sin{\theta}+\sigma_{r}\cos{\theta}-\sigma_{\theta}\cos{\theta}-2\tau_{\theta r}\sin{\theta}\right)+\cos{\theta}\frac{\partial\tau_{rz}}{\partial z}-\frac{\partial\tau_{\theta z}}{\partial z}\sin{\theta} \end{eqnarray}

同様に、$y, z$軸方向の運動方程式である式(2), (3)についても式(10), (13), (15)を代入して整理します。式(2)については

\begin{eqnarray} \rho\frac{\partial^2u_r}{\partial t^2}\sin{\theta}+\rho\frac{\partial^2u_\theta}{\partial t^2}\cos{\theta}&=& \left(\cos{\theta}\frac{\partial}{\partial r}-\frac{\sin{\theta}}{r}\frac{\partial}{\partial\theta}\right)\left\{\left(\sigma_{r}-\sigma_{\theta}\right)\sin{\theta}\cos{\theta}+\tau_{\theta r}\cos{2\theta}\right\}+\left(\sin{\theta}\frac{\partial}{\partial r}+\frac{\cos{\theta}}{r}\frac{\partial}{\partial\theta}\right)\left(\sigma_{r}\sin^2{\theta}+\tau_{\theta r}\sin{2\theta}+\sigma_{\theta}\cos^2{\theta}\right)+\frac{\partial}{\partial z}\left(\tau_{rz}\sin{\theta}+\tau_{\theta z}\cos{\theta}\right) \nonumber\\ \rho\frac{\partial^2u_r}{\partial t^2}\sin{\theta}+\rho\frac{\partial^2u_\theta}{\partial t^2}\cos{\theta}&=& \frac{\partial\sigma_{r}}{\partial r}\sin{\theta}+\frac{\partial\tau_{\theta r}}{\partial r}\cos{\theta}+\frac{1}{r}\left(\frac{\partial\tau_{\theta r}}{\partial\theta}\sin{\theta}+\frac{\partial\sigma_{\theta}}{\partial\theta}\cos{\theta}+\sigma_{r}\sin{\theta}-\sigma_{\theta}\sin{\theta}+2\tau_{\theta r}\cos{\theta}\right)+\frac{\partial\tau_{rz}}{\partial z}\sin{\theta}+\frac{\partial\tau_{\theta z}}{\partial z}\cos{\theta} \end{eqnarray}

となり、式(3)については

\begin{eqnarray} \rho\frac{\partial^2w}{\partial t^2}&=& \left(\cos{\theta}\frac{\partial}{\partial r}-\frac{\sin{\theta}}{r}\frac{\partial}{\partial\theta}\right)\left(\tau_{rz}\cos{\theta}-\tau_{\theta z}\sin{\theta}\right)+\left(\sin{\theta}\frac{\partial}{\partial r}+\frac{\cos{\theta}}{r}\frac{\partial}{\partial\theta}\right)\left(\tau_{rz}\sin{\theta}+\tau_{\theta z}\cos{\theta}\right)+\frac{\partial\sigma_z}{\partial z} \nonumber\\ \rho\frac{\partial^2w}{\partial t^2}&=&\frac{1}{r}\frac{\partial}{\partial r}\left(r\tau_{rz}\right)+\frac{1}{r}\frac{\partial\tau_{\theta z}}{\partial\theta}+\frac{\partial\sigma_z}{\partial z} \end{eqnarray}

となります。しかしここで、最終的に求めたい運動方程式は次のように円筒座標系での各軸方向の式であるため3

\begin{equation} \rho\frac{\partial^2u_r}{\partial t^2}=\cdots ,\ \ \rho\frac{\partial^2u_\theta}{\partial t^2}=\cdots,\ \ \rho\frac{\partial^2w}{\partial t^2}= \cdots \nonumber \end{equation}

式(16), (17)では不十分なことが分かります。そこで、いったん式(16)+式(17)とします。

\begin{equation} \frac{\partial^2u_r}{\partial t^2}\left(\cos{\theta}+\sin{\theta}\right)+\rho\frac{\partial^2u_\theta}{\partial t^2}\left(\cos{\theta}-\sin{\theta}\right)=\frac{\partial\sigma_{r}}{\partial r}\left(\cos{\theta}+\sin{\theta}\right)+\frac{\partial\tau_{\theta r}}{\partial r}\left(\cos{\theta}-\sin{\theta}\right)+\frac{1}{r}\left\{\frac{\partial\tau_{\theta r}}{\partial\theta}\left(\cos{\theta}+\sin{\theta}\right)+\frac{\partial\sigma_{\theta}}{\partial\theta}\left(\cos{\theta}-\sin{\theta}\right)+\sigma_{r}\left(\cos{\theta}+\sin{\theta}\right)-\sigma_{\theta}\left(\cos{\theta}+\sin{\theta}\right)+2\tau_{\theta r}\left(\cos{\theta}-\sin{\theta}\right)\right\}+\frac{\partial\tau_{rz}}{\partial z}\left(\cos{\theta}+\sin{\theta}\right)+\frac{\partial\tau_{\theta z}}{\partial z}{\left(\cos{\theta}-\sin{\theta}\right)} \end{equation}

ここで分かる通り、式(19)は$\cos{\theta}+\sin{\theta}$の項と$\cos{\theta}-\sin{\theta}$の項に分けられます。そして、次の関係を用い式(19)を変形します。

\begin{equation} \cos{\theta}+\sin{\theta}=\frac{2}{\sqrt2}\left(\sin{\frac{\pi}{4}}\cdot\cos{\theta}+\sin{\theta}\cdot\cos{\frac{\pi}{4}}\right)=\frac{2}{\sqrt2}\sin{\left(\frac{\pi}{4}+\theta\right)}\\ \cos{\theta}-\sin{\theta}=\frac{2}{\sqrt2}\left(\sin{\frac{\pi}{4}}\cdot\cos{\theta}-\sin{\theta}\cdot\cos{\frac{\pi}{4}}\right)=\frac{2}{\sqrt2}\sin{\left(\frac{\pi}{4}-\theta\right)} \nonumber \end{equation}

これにより、

\begin{equation} \frac{\partial^2u_r}{\partial t^2}\sin{\left(\frac{\pi}{4}+\theta\right)}+\rho\frac{\partial^2u_\theta}{\partial t^2}\sin{\left(\frac{\pi}{4}-\theta\right)}=\frac{\partial\sigma_{r}}{\partial r}\sin{\left(\frac{\pi}{4}+\theta\right)}+\frac{\partial\tau_{\theta r}}{\partial r}\sin{\left(\frac{\pi}{4}-\theta\right)}+\frac{1}{r}\left\{\frac{\partial\tau_{\theta r}}{\partial\theta}\sin{\left(\frac{\pi}{4}+\theta\right)}+\frac{\partial\sigma_{\theta}}{\partial\theta}\sin{\left(\frac{\pi}{4}-\theta\right)}+\sigma_{r}\sin{\left(\frac{\pi}{4}+\theta\right)}-\sigma_{\theta}\sin{\left(\frac{\pi}{4}+\theta\right)}+2\tau_{\theta r}\sin{\left(\frac{\pi}{4}-\theta\right)}\right\}+\frac{\partial\tau_{rz}}{\partial z}\sin{\left(\frac{\pi}{4}+\theta\right)}+\frac{\partial\tau_{\theta z}}{\partial z}{\sin{\left(\frac{\pi}{4}-\theta\right)}} \end{equation}

となります。上式の両辺で$\sin{\left(\frac{\pi}{4}+\theta\right)}$と$\sin{\left(\frac{\pi}{4}-\theta\right)}$の項の係数は等しいはずなので、

\begin{equation} \rho\frac{\partial^2u_r}{\partial t^2}=\frac{\partial\sigma_{r}}{\partial r}+\frac{1}{r}\left(\frac{\partial\tau_{\theta r}}{\partial\theta}+\sigma_{r}-\sigma_{\theta}\right)+\frac{\partial\tau_{rz}}{\partial z}\\ \rho\frac{\partial^2u_\theta}{\partial t^2}=\frac{\partial\tau_{\theta r}}{\partial r}+\frac{1}{r}\left(\frac{\partial\sigma_{\theta}}{\partial\theta}+2\tau_{\theta r}\right)+\frac{\partial\tau_{\theta z}}{\partial z} \end{equation}

が成り立ちます。最終的に、円柱座標系での運動方程式を以下の通り表せます。

\begin{equation} \left\{ \begin{array}{c} \rho\frac{\partial^2u_r}{\partial t^2}=\frac{\partial\sigma_{r}}{\partial r}+\frac{1}{r}\left(\frac{\partial\tau_{\theta r}}{\partial\theta}+\sigma_{r}-\sigma_{\theta}\right)+\frac{\partial\tau_{rz}}{\partial z} \\ \rho\frac{\partial^2u_\theta}{\partial t^2}=\frac{\partial\tau_{\theta r}}{\partial r}+\frac{1}{r}\left(\frac{\partial\sigma_{\theta}}{\partial\theta}+2\tau_{\theta r}\right)+\frac{\partial\tau_{\theta z}}{\partial z}\\ \rho\frac{\partial^2w}{\partial t^2}=\frac{1}{r}\frac{\partial}{\partial r}\left(r\tau_{rz}\right)+\frac{1}{r}\frac{\partial\tau_{\theta z}}{\partial\theta}+\frac{\partial\sigma_z}{\partial z} \end{array} \right. \end{equation}

最後に

特に式(19)以降の式展開については解説している資料が見つからなかったので自分なりにまとめてみました。ちなみにですが、慣性力を体積力に書き換えることで静解析に必要な平衡方程式を導くことも可能です


  1. 荻博次、”弾性力学”、共立出版 (2011). ↩︎

  2. 詳しい話はこちらで説明しています。 ↩︎

  3. 今回の記事では説明しませんが、最終的に振動数方程式を導く際に各軸方向の運動方程式を導いている必要があります。 ↩︎