モチベーション

流体力学を学ぶ者として、境界層方程式は使えないとでしょ(素人

最終アウトプット

境界層内の速度分布(この記事)
境界層内の温度分布(今後)
マッハ数による変化(今後)

1. 境界層方程式とは？

境界層方程式を語る前にまずはナビエストークス方程式を軽くおさらい。
ナビエストークス方程式は流体の運動を記述するスーパーな方程式だが、取り扱いに大きな制約がある。
非圧縮性NS式（Navier–Stokes equations）
$$ \begin{align}
\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\nu \frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\frac{\mu}{\rho}(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2})\
\end{align}$$

$$ \begin{align}
\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial x}+\nu \frac{\partial v}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\frac{\mu}{\rho}(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}) \
\end{align}$$

この方程式は、$x$と$y$の2階偏微分がある、いわゆる楕円型の偏微分方程式である。
方程式を解くには、初期条件と$x$(流れ方向)と$y$にそれぞれ2つの境界条件が必要となる。
（すなわち上流と下流の影響を受ける）
実際に解く場合は、空間・時間を離散化する必要があるが、精度を求めて空間/空間を細かく分割すると計算コストが膨大になり、昔は解くことが難しかった。

そこで登場したのが境界層方程式である。
境界層方程式は、壁近傍の流れだけに粘性を考慮することでNS式を簡略化している。
式はこんな感じ。

\begin{cases}
\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\nu\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\frac{\mu}{\rho}\frac{\partial^2 u}{\partial y^2}\\
\frac{\partial p}{\partial y} =0\\
\frac{\partial u}{\partial x}+\frac{\partial \nu}{\partial y} =0
\end{cases}

この方程式の2階微分は$y$の2階偏微分のみで、いわゆる放物型の偏微分方程式である。
この偏微分方程式を解くためには、初期条件と$x$に対する1つの境界条件、$y$に対する2つの境界条件があればよい。（すなわち、上流境界条件のみで解が求められ、下流方向へは積分を進めることができる。

2. 境界層方程式の利用

今は計算機の性能も上がり、100年前ほど制約は厳しくない。
「流れ場の情報が欲しければ、そこらへんに転がっている(※環境による)RANSとかで解けばいいじゃん、なぜ境界層方程式を使うの？」という意見も当然だ。
それでも境界層方程式は、境界層内の速度、濃度や温度の分布をパパっと求めることができるため、今でも様々なアプリケーションに使われている。
最近も風洞模型の検討で境界層内の速度分布が必要だったので使った。
備忘録もかねて、境界層方程式を使って平板境界層内の速度分布を求めてみる。
ここでは、境界層方程式の詳細な導出までは行わないので悪しからず。

前提条件

平板と並行方向に$x$軸方向、垂直方向に$y$軸をとる
流れは平板と並行の向きに流れ、逆流、再循環および剥離がないと仮定する
流れ場は十分発達し、定常であると仮定する
$x$方向速度$u$は、平板の長さ方向($x$方向)に以下の相似形の速度分布を持つと仮定する。
$$\frac{u}{U_\infty}=\phi_1 (\frac{y}{\delta}) \tag{2-1}$$
境界層の厚みは、以下の特徴を持つ
- 平板端からの距離$x$と動粘度$\nu$の積の平方根に比例する
- 主流速度$U_\infty$の平方根に反比例する

ηと流れ関数ψの導入

境界層の厚みについての仮定より、
$$\frac{u}{U_\infty}=\phi_2 (y\sqrt{\frac{U_\infty}{\nu x}}) \tag{2-2}$$
と書ける。
次に、相似変数$\eta$を導入する。

\eta\equiv y\sqrt{\frac{U_\infty}{\nu x}} \tag{2-3}

相似解を以下と仮定する。

\frac{u}{U} = g(\eta)\tag{2-4}

次に、流れ関数 $\Psi$ を導入すると、各方向速度との関係は以下。

u\equiv\frac{\partial \Psi}{\partial y} \tag{2-5} \\

v\equiv-\frac{\partial \Psi}{\partial x}\tag{2-6}

なお、$\Psi$は自動的に連続の式を満足する。←これが流れ関数を導入した最大のメリット

境界層方程式の無次元化

境界層方程式（x方向NS式, y方向NS, 連続の式)

\begin{cases}
\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\nu\frac{\partial u}{\partial y} = -\frac{1}{\rho}\frac{\partial p}{\partial x}+\frac{\mu}{\rho}\frac{\partial^2 u}{\partial y^2}\\
\frac{\partial p}{\partial y} =0\\
\frac{\partial u}{\partial x}+\frac{\partial \nu}{\partial y} =0 \tag{2-7} 
\end{cases}

式(2-5)を$y$で積分すると、

\begin{align}
\Psi &= \int udy + c(y) \\
 &=\int U_\infty g(\eta) dy + c(y) \\
 &=U_\infty \int g(\eta)\frac{dy}{d\eta}d\eta + c(y)
\end{align}

式(2-3)より、$\frac{dy}{d\eta} = \sqrt{vx/U_0}$となるため、

\begin{align}
\Psi &= U_\infty \sqrt{vx/U_\infty} \int g(\eta)d\eta + c(y) \\ 
 &= \sqrt{v x  U_\infty} \int g(\eta) d\eta + c(y) \tag{2-8}
\end{align}

式(2-6)に式(2-8)を代入すると、

v = -\frac{\partial(\sqrt{v x U_\infty} \int g(/eta)d\eta)}{\partial x}+\frac{\partial(c(y))}{\partial x}

仮定より、$y=0$ で$v=0$の条件を満たすためには、$c(y)=0$とすればよい。
したがって、$\int g(\eta)d\eta = f(\eta)$とおくと、
流れ関数$\Psi$は
$$\Psi=\sqrt{\nu x U_\infty}f(\eta) \tag{2-9}$$
となる。

式(2-5)より、

\begin{align}
u &=\frac{\partial}{\partial y}\sqrt{\nu x U_\infty}f(\eta) \\ 
&= \sqrt{\nu x U_\infty}\frac{\partial f}{\partial \eta}\frac{\partial \eta}{\partial y}\\  
&= \sqrt{\nu x U_\infty} \sqrt{\frac{U_\infty}{\nu x}}\frac{\partial f}{\partial \eta}\\
&= U_\infty f' \tag{2-10}
\end{align}

式(2-6)より、

v = -1/2 \sqrt{\frac{vU_\infty}{x}}f(\eta)-\sqrt{vU_\infty x}f'\frac{\partial \eta}{\partial x} \\

ここで、式(2-3)より、

\begin{align}
\frac{\partial \eta}{\partial x} &= -\frac{\eta}{2 x} \\
\frac{\partial \eta}{\partial y} &= \frac{1}{\sqrt{vx/U_\infty}} \tag{2-11}
\end{align}

なので、

\begin{align}
v &= -1/2 \sqrt{\frac{vU_\infty}{x}}f(\eta)+\sqrt{vU_\infty x}f'\frac{\eta}{2x} \\
 &= \frac{1}{2}\sqrt{\frac{\nu U_\infty}{x}}(\eta f'-f) \tag{2-12}
\end{align}

よって式(2-10)より、

\begin{align}
\frac{\partial u}{\partial x} &= \frac{\partial u}{\partial \eta} \frac{\partial \eta}{\partial x} = -\frac{\eta U_\infty}{2x} f'' \\ 

\frac{\partial u}{\partial y} &= \frac{\partial u}{\partial \eta} \frac{\partial \eta}{\partial y} = U_\infty f'' \sqrt{\frac{U_\infty}{vx}} \\ 

\frac{\partial^2 u}{\partial y^2} &= \frac{\partial}{\partial \eta} (U_\infty f'' \sqrt{\frac{U_\infty}{vx}} \frac{\partial \eta}{\partial y} 
= \frac{U_\infty^2}{\nu x}f'''   \tag{2-13}
\end{align}

となる。

連立常微分方程式の定式化

式(2-7) を変形する。
x方向に流れが変化しない、かつ定常の流れという仮定より、

\begin{align}
\frac{\partial u}{\partial t} &= 0 \\
\frac{\partial p}{\partial x} &= 0 \tag{2-14}
\end{align}

となり、最終的には以下。
$$2f''' + f f'' = 0 \tag{2-15}$$

対象とする境界層では、

壁面($y = 0$)では各速度成分が0
壁面から十分離れた位置($y=\infty$)では、x方向流速が主流流速と等しい。
境界条件として整理すると以下の通り。

\begin{cases}
u(x,0) = v(x,0) = 0 \\
u(x,\infty) = U_\infty
\end{cases}

$y = 0$のとき、$\eta = 0$
$y = \infty$のとき、$\eta = \infty$
境界条件を書き換えると、

\begin{cases}
\eta=0 &:f&=f'=0 \\ 
\eta=\infty &:f'&=1 \\ 
\end{cases}

以上で境界層方程式を解くための下準備は完了。
やっと定式化に入る。
境界層方程式を以下の連立常微分方程式に書き換える。（いずれも$\eta$の関数）

\begin{cases}
y_0' &= y_1 \\ 
y_1' &= y_2 \\ 
y_2' &= -\frac{1}{2}y_0 y_2
\end{cases}

境界条件も書き換える

\begin{cases}
y_0(0) &= 0 \\ 
y_1(0) &= 0 \\ 
y_2(\infty) &=1 \\ 
\end{cases}

これで境界層方程式の定式化が完了。
この問題は、$\eta=0$で$y_2$ の値が与えられておらず、代わりに$\eta=\infty$で$y_1=1$が与えられている。
この問題は境界値問題として解く必要がある。
境界値問題を解くには、scipyのslove_bvpが使える。

3. Pythonでの境界層方程式の求解

必要なモジュールをインポートする

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
from sympy import *

境界層方程式の定義

$y_0~y_2$を配列yに格納
$\eta$は変数
$y_0~y_2$は、ひとつ前の値を用いて計算する

def fun(η, y):
    dy0dη = y[1]
    dy1dη = y[2]
    dy2dη = -0.5 * y[0] * y[2]
    return np.array([dy0dη, dy1dη, dy2dη])

境界条件の定義

yaの変数の定義域の片側を表す。-> $y_0,y_1=0$@$\eta=0$
ybは変数の定義域のもう片側を表す -> $y_1=0$@$\eta=0$
返り値は、各境界値の値との差分を返す。

def bc(ya, yb):
    return np.array([ya[0]-0, ya[1]-0, yb[1] - 1])

変数 𝜂の定義

変数 𝜂の定義域は、境界条件の収束判定ができるくらいには広くとる必要がある。
今回は、 𝜂=10 とした。
次に境界値問題ではあるけれど、初期値の設定する。初期値と言っても各 𝜂
での方程式の値が格納する変数である。

# 初期値
η_values = np.linspace(0, 10, 101)
y_guess = np.zeros((3, η_values.size))

連立常微分方程式を解きます。

やっと来ました、solve_bvpを用い、微分方程式を解いてみます。

# solve_bvpを使用して連立常微分方程式を解く
solution = solve_bvp(fun, bc, η_values, y_guess)
# 結果評価用のy1の厳密解
Howarth= np.loadtxt("Howarth.csv", delimiter=',')
#
# 解をプロット
plt.plot(solution.x, solution.y[0], label='y0(η)')
plt.plot(solution.x, solution.y[1], label='y1(η)')
plt.plot(solution.x, solution.y[2], label='y2(η)')
plt.scatter(Howarth[:,0:1], Howarth[:,1:2], label='Howarth')

plt.xlabel('η')
plt.ylim(0,1.2)
plt.xlim(0,10)
plt.legend()
plt.show()

プログラムの出力はコチラ！
𝑦1 の出力が厳密解(Howwarth)と一致していることがわかる。

では、問題を設定したときにはわからなかった、 𝜂=0の時の 𝑦2 を確認してみると、
$$𝑦2 = 0.33205$$
参考文献でShooting Methodで計算した結果は、$𝑦2 = 0.33206$ なので概ね正しいね。

まとめ

この記事では、
1 境界層方程式の紹介
1 微分方程式としての整理
1 pythonを使った微分方程式の求解
を紹介しました。次回は温度分布を求めてみたいとおもいます。

参考文献

https://www.jstage.jst.go.jp/article/htsj1999/43/178/43_178_9/_pdf