ナビエ-ストークス方程式は一般には解析的に解くことができません（手計算で解けない）。しかし、流れが十分に単純な場合には解析解を求めることができ、流体の運動を理解する上で非常に役立つ知見を得られます。なぜなら、現実の複雑な粘性流も局所的には解析解が存在するような基本的な流れに帰着できることが多いからです。

この記事では、非圧縮性ナビエ-ストークス方程式から出発し、3次元流れであるHagen-Poiseuille流の解析解を導きます。元は偏微分方程式ですが、連鎖律、偏微分の理論を駆使してシンプルな常微分方程式を導出し、微分方程式の解法で一般解を求めていきます。大学の前半で学んだ知識をフル動員し、工学的に重要な法則を導きます。

ナビエ-ストークス方程式

非圧縮性粘性流体に対するナビエ-ストークス方程式は次の形で表されます。

$$\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \cdot \nabla \boldsymbol{u} = -\frac{1}{\rho} \nabla p + \nu \Delta \boldsymbol{u} + \boldsymbol{f} \tag{1}$$

平行流

流線がすべて定直線に平行な直線である流れを平行流といい、Hagen-Poiseuille流は平行流の一種です。定直線の方向に直交座標$(x, y, z)$の$x$軸を取ると、平行流の定義からその速度は$\boldsymbol{u} = (u, 0, 0)$で表されます。$v=w=0$、非圧縮性流体の連続の式より

$$\frac{\partial u}{\partial x} = 0 \tag{2}$$

つまり、$u$は$x$軸方向に一定です。

また、外力がないものとすれば、$v=w=0$を式$(1)$の$y$、$z$成分に代入することで

$$\frac{\partial p}{\partial y} = \frac{\partial p}{\partial z} = 0 \tag{3}$$

となり、$p$は$y$軸、$z$軸方向に一定です。

さらに、式$(1)$の$x$成分は式$(2)$より

$$\frac{\partial u}{\partial t} = -\frac{1}{\rho} \frac{\partial p}{\partial x} + \nu(\frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}) \tag{4}$$

となります。

$u$を含む項は$x$によらず、$p$を含む項は$y$、$z$によらないので、これらはそれぞれ時間$t$だけで表すことができるはずです。

$$\frac{\partial u}{\partial t} - \nu(\frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2}) = -\frac{1}{\rho} \frac{\partial p}{\partial x} = \alpha(t)\tag{5}$$

定常流なら、

$$\frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} = \frac{1}{\mu} \frac{\partial p}{\partial x} = -\frac{\alpha}{\nu} \tag{6}$$

式$(6)$をベースにHagen-Poiseuille流の解析解を導出していきます。

Hagen-Poiseuille流

円管断面の管を通る流れを考えます。平行流の定直線方向に$x$軸をとれば、$y$軸、$z$軸は下図の左のようになります。この流れの流速分布は同心円状になると期待できるので、直交座標よりも極座標を使った方が見通しがよさそうです。そこで、下図の右のように$(y, z)$を極座標$(r, \theta)$へと変換します。

式$(6)$の左辺を変換するためには、$\partial^2 u/\partial y^2$、$\partial^2 u/\partial z^2$を$r$と$\theta$の微分で表現する必要があります。そのために、連鎖律を利用します。

いま、流速$u$は$r$と$\theta$で表すことができるはずであり、$r(y, z)$、$\theta(y, z)$であるので、$u$は$y$、$z$の合成関数と見立てることができます。連鎖律に当てはめると

$$\begin{align} \frac{\partial u}{\partial y} &= \frac{\partial u}{\partial r} \frac{\partial r}{\partial y} + \frac{\partial u}{\partial \theta} \frac{\partial \theta}{\partial y} = \frac{\partial u}{\partial r} \frac{\partial r}{\partial y} \notag \\ \frac{\partial u}{\partial z} &= \frac{\partial u}{\partial r} \frac{\partial r}{\partial z} + \frac{\partial u}{\partial \theta} \frac{\partial \theta}{\partial z} = \frac{\partial u}{\partial r} \frac{\partial r}{\partial z} \tag{7} \end{align}$$

2回目の式変形には、流速分布が同心円状になる、すなわち$\partial u / \partial \theta = 0$とする仮定を利用しました。

続いて、式$(7)$の$\partial r/\partial y$と$\partial r/\partial z$は以下のように$\theta$のみの関数に式変形できます。

$$\begin{align} \frac{\partial r}{\partial y} &= \frac{\partial}{\partial y} \left\{ (y^2 + z^2)^{1/2}\right\} = \frac{y}{r} = \cos \theta \notag \\ \frac{\partial r}{\partial z} &= \frac{\partial}{\partial z} \left\{ (y^2 + z^2)^{1/2}\right\} = \frac{z}{r} = \sin \theta \tag{8} \end{align}$$

よって式$(7)$は以下のように$r$と$\theta$のみで表すことができます。

$$\begin{align} \frac{\partial u}{\partial y} &= \frac{\partial u}{\partial r} \cos \theta \notag \\ \frac{\partial u}{\partial z} &= \frac{\partial u}{\partial r} \sin \theta \tag{9} \end{align}$$

目標は2階微分を$r$と$\theta$のみで表すことなので、同じ要領で以下のとおり計算を進めます。

$$\begin{align} \frac{\partial^2 u}{\partial y^2} &=\frac{\partial}{\partial y}\left(\frac{\partial u}{\partial y}\right) \notag\\ &= \frac{\partial}{\partial r}\left(\frac{\partial u}{\partial y}\right)\frac{\partial r}{\partial y} + \frac{\partial}{\partial \theta}\left(\frac{\partial u}{\partial y}\right)\frac{\partial \theta}{\partial y} \notag\\ &= \frac{\partial}{\partial r} \left( \frac{\partial u}{\partial r}\cos\theta \right) \frac{\partial r}{\partial y} + \frac{\partial}{\partial \theta} \left( \frac{\partial u}{\partial r}\cos\theta \right) \frac{\partial \theta}{\partial y} \notag\\ &= \left( \frac{\partial ^2 u}{\partial r^2}\cos\theta \right)\cos\theta + \left( -\frac{\partial u}{\partial r}\sin\theta \right) \left( -\frac{\sin\theta}{r} \right) \notag\\ &= \frac{\partial^2 u}{\partial r^2}\cos^2\theta + \frac{\sin^2\theta}{r}\frac{\partial u}{\partial r} \tag{10} \end{align}$$

同様に、$z$についての2階微分も計算します。

$$\begin{align} \frac{\partial^2 u}{\partial z^2} &=\frac{\partial}{\partial z}\left(\frac{\partial u}{\partial z}\right) \notag\\ &= \frac{\partial}{\partial r}\left(\frac{\partial u}{\partial z}\right)\frac{\partial r}{\partial z} + \frac{\partial}{\partial \theta} \left( \frac{\partial u}{\partial z} \right) \frac{\partial \theta}{\partial z} \notag\\ &= \frac{\partial}{\partial r} \left( \frac{\partial u}{\partial r}\sin\theta \right) \frac{\partial r}{\partial z} + \frac{\partial}{\partial \theta} \left( \frac{\partial u}{\partial r}\sin\theta \right) \frac{\partial \theta}{\partial z} \notag\\ &= \left( \frac{\partial ^2 u}{\partial r^2}\sin\theta \right)\sin\theta + \left( \frac{\partial u}{\partial r}\cos\theta \right) \left( \frac{\cos\theta}{r} \right) \notag\\ &= \frac{\partial^2 u}{\partial r^2}\sin^2\theta + \frac{\cos^2\theta}{r}\frac{\partial u}{\partial r} \tag{11} \end{align}$$

式$(10)$と式$(11)$を用いることで、式$(6)$は次式のように極座標へ変換することができます。2変数$(y, z)$から1変数$r$の常微分方程式へ帰着できています。 このおかげで手計算で解くことができます。

$$\frac{d^2 u}{d r^2} + \frac{1}{r}\frac{d u}{d r} = - \frac{\alpha}{\nu} \tag{12}$$

ここでポイントなのは、式$(12)$の左辺が積の微分を展開したような形になっていると気付くことです。$v = du/dr$とおいてみると式$(12)の左辺は

$$\frac{dv}{dr} + \frac{1}{r}v \tag{13}$$

となります。この形から、$rv$を積の微分で展開したものと予想して試してみます。

$$\frac{d}{dr}(rv) = r\frac{dv}{dr} + v = r \left(\frac{dv}{dr} + \frac{1}{r} v \right) \tag{14}$$

式$(13)$と式$(14)$を見比べれば、式$(12)$の右辺は

$$\frac{d^2 u}{d r^2} + \frac{1}{r}\frac{d u}{d r} = \frac{dv}{dr} + \frac{1}{r}v = \frac{1}{r}\frac{d}{dr}\left(r \frac{du}{dr} \right) \tag{15}$$

と変形できることに気付けるはずです。式$(12)$を眺めて、積の微分の形っぽいなと感じることが肝要です。

式$(12)$と式$(15)$より、

$$\begin{align} \frac{d}{dr}\left(r \frac{du}{dr} \right) &= -\frac{\alpha}{\nu}r \notag \\ r \frac{du}{dr} &= -\frac{\alpha}{2\nu}r^2 + C_1 \notag \\ \frac{du}{dr} &= -\frac{\alpha}{2\nu}r + \frac{C_1}{r} \tag{16} \end{align}$$

さらに積分すると

$$u = -\frac{\alpha}{4\nu}r^2 + C_1 \ln r + C_2 \tag{17}$$

が得られます。ここで、$r=0$（管の中心）で流速$u$が有限であるためには、$C_1 = 0$でなければなりません。$C_1 \neq 0$だと$r \to 0$で$\ln r \to -\infty$となり、流速が負の無限大となってしまい、物理的に不合理だからです。よって

$$u = -\frac{\alpha}{4\nu}r^2 + C_2 \tag{18}$$

が得られます。いまは円管を考えているので、内径を$a$として管壁（$r=a$）で$u=0$（滑りなし条件）を適用すると

$$0 = -\frac{\alpha}{4\nu}a^2 + C_2 \quad \therefore C_2 = \frac{\alpha}{4\nu}a^2 \tag{19}$$

したがって、速度分布$u$は

$$u = \frac{\alpha}{4\nu}(a^2 - r^2) \tag{20}$$

となります。

これがHagen-Poiseuille流の速度分布です。2次元Poiseuille流と同様に放物線状の分布となっています。

最大流速と流量

式$(20)$から明らかなように、流速は$r = 0$（管の中心）で最大となり、その値は

$$u_{\max} = \frac{\alpha a^2}{4\nu} = -\frac{a^2}{4\mu}\frac{\partial p}{\partial x} \tag{21}$$

となります。また、体積流量$Q$は

$$Q = \int_0^a u \cdot 2\pi r \, dr = \frac{\pi \alpha a^4}{8\nu} \tag{22}$$

つまり、ある配管を流体が流れるとき、その流量は半径の4乗に比例し、圧力勾配に比例します。これがHagen-Poiseuilleの法則です。半径が2倍になれば流量は16倍になるという、工学的に非常に重要な法則です。

Hagen-Poiseuille流は粘性流の式$(6)$の厳密解なので、すべての$Re$で成り立つはずです。しかし、平均流速がある限界を超え、$Re$が2000を超えたあたりで、流れは乱流に移行します。この話はまたの機会に。

参考図書

• 巽友正「流体力学」（培風館）：より厳密な数学的定式化と幅広いトピックをカバーしており、大学院レベルの学習に適しています。読むには微積分、ベクトル解析の大学レベルの知識が必要です。