二次元流れ場


二次元において、非圧縮性流れのナビエ・ストークス方程式
$$
\frac{\partial \vec{u}}{\partial t} + (\vec{u}\cdot\nabla)\vec{u}= \frac{1}{\rho}\nabla p + \frac{\mu}{\rho} \nabla^2 \vec{u} + \vec{K}
$$
流体にかかる重力を\(K=k\hat{x}\)とする。
圧力場を\(p(x,y)=p_0/\sqrt{x^2+y^2}\)とし、定常状態を考える。
$$
(\vec{u}\cdot\nabla)\vec{u}= \frac{1}{\rho}\nabla p + \frac{\mu}{\rho} \nabla^2 \vec{u} + \vec{K}
$$
$$
(u\frac{\partial}{\partial x}+v\frac{\partial}{\partial y})\vec{u}= \frac{1}{\rho}\nabla p + \frac{\mu}{\rho} \nabla^2 \vec{u} + k\hat{x}
$$

(2次元)直管損失


負荷のかかる流体計算を簡略化するための常套手段に等価回路近似がある。まず、定常状態において粘性抵抗を求め、非定常状態においてはどの程度その近似が成り立っているのかを検証する。

前回の記事において(無限級数の形であるが)解析的にポワズイユ流の非定常解を計算した。今回は、前回と同様の以下の流れにおいてy方向の分布は考えず、1次元的な流れの式に変えて考える。
poiseuille
もともとの微分方程式
$$
\frac{\partial u}{\partial t} = \frac{\Delta p}{\rho l} + \frac{\mu}{\rho} \frac{\partial^2 u}{\partial y^2}
$$
を\(y\)方向の平均をとるため積分すると、
$$
\frac{\partial \bar{u}}{\partial t} = \frac{\Delta p}{\rho l} + \frac{\mu}{\rho} \frac{1}{d}\left[ \frac{\partial u}{\partial y} \right]_{y=0}^d
$$
と表される。ここで、右辺第一項により流体は圧力による加速を受け、右辺第二項により粘性拡散される。各項の働きを考えれば左辺はコイル、右辺第一項は印加電圧、右辺第二項を抵抗というように置き換えて考えることができそうである。抵抗は電流に比例して働くのと同様、右辺第二項が\(u\)に比例するとして式を書き直したい。
本来であれば、右辺第二項を、非定常の解析解を用いて書き直したいところだが、無限級数の形となっているためそれができない。今回は、 前々回求めた定常解をもとに考えてみる。
定常解\(u(y)=-\frac{\Delta p}{2\mu l}y(y-d)\)を用いると右辺第二項は
$$
\frac{\mu}{\rho} \frac{1}{d}\left[ \frac{\partial u}{\partial y} \right]_{y=0}^d=\frac{\mu}{\rho}\frac{1}{d}\left(-\frac{\Delta p}{2\mu l}\right)\left[ 2y-d \right]_{y=0}^d=-\frac{\mu}{\rho}\frac{\Delta p}{\mu l}
$$
となる。右辺第二項を\(\bar{u}\)に比例する形にするため、定常解の\(y\)方向平均をとり
$$
\bar{u}=-\int_0^d\frac{\Delta p}{2\mu l}y(y-d) dy = \frac{1}{12}\frac{\Delta p}{\mu l}d^2
$$
を得る。これより\(\frac{\Delta p}{\mu l}=\frac{12}{d^2}\bar{u}\)の関係を適用させると、\(t\)だけの全微分方程式
$$
\frac{\partial \bar{u}_{\rm app}}{\partial t} = \frac{\Delta p}{\rho l} – \frac{\mu}{\rho} \frac{12}{d^2}\bar{u}_{\rm app}
$$
が完成した。あくまでもこの式は非定常状態を表すために粘性抵抗を定義した。以下、これが非定常状態においてはどの程度正しい近似になっているのかを考える。

検証

時間発展

この微分方程式の時間発展がどの程度もとの式と違いがあるのかを検証するため、前回の記事と同様、\(t=0\)で\(\bar{u}_{\rm app}=0\)とした初期条件で近似解を得てみる。右辺全体を\(X\)とおくと
$$
\frac{dX}{dt} =-\frac{12\mu}{\rho d^2} X
$$
より、\(X=C\exp\left(-\frac{12\mu}{\rho d^2} t\right)\)となる。
初期条件も考えると\(C\)も決定され解は
$$
\bar{u}_{\rm app}=\frac{\Delta p d^2}{12\mu l}\left( 1 – \exp \left(- \frac{12\mu}{\rho d^2} t \right) \right)
$$
となる。

解析解との比較

gnuplotを用いて、数値的に解析解
$$
\bar{u}(t)=\sum_{n=1}^\infty \frac{8 d^2\Delta p}{\mu l(2n-1)^4\pi^4} \left( 1- \exp\left( -\frac{\mu (2n-1)^2\pi^2}{\rho d^2} t \right) \right)
$$
と比較を行うと、以下の図の通りになっており、良い近似になっていることがわかる。

comparison

同様の図はgnuplotにて以下のように入力することで比較可能だ。

補正

過渡現象開始から収束までの全体を合わせる補正

この図を見ると、時定数を補正すればさらに良い近似になりそうに思われる。そこで、補正係数\(C_{\rm cor}\)を定義して
$$\frac{\partial \bar{u}_{\rm app}}{\partial t} = C_{\rm cor} \left( \frac{\Delta p}{\rho l} – \frac{\mu}{\rho} \frac{12}{d^2}\bar{u}_{\rm app} \right)$$
としたときの解
$$\bar{u}_{\rm app}=\frac{\Delta p d^2}{12\mu l}\left( 1 – \exp \left(- C_{\rm cor} \frac{12\mu}{\rho d^2} t \right) \right)$$
を解析解にフィッティングした。その結果、補正係数は\(C_{\rm cor}=0.834(1)\)という結果になり、比較のため再度プロットすると下図のようになった。
approximation-with-correctionかなりよく一致していることが確認できる。丸め誤差等もあるのでどこまでgnuplotで詳細に評価できるのか不明であるが、近似解から解析解を引いた後、解析解で規格化してプロットした結果は以下の図となった。
appro-error過渡現象における\(t=0\)の近傍においては誤差が大きそうである。誤差の大きさを見積もると、

過渡現象直近に対する補正

そこで、\(t=0\)における\(\bar{u}\)を合わせるように補正をし直してみる。\(t\)の1次までテーラー展開すると、解析解は
$$
\begin{eqnarray}
\bar{u}(t)&\sim&\sum_{n=1}^\infty \frac{8 \Delta p}{\rho l(2n-1)^2\pi^2} t = \frac{\Delta p}{\rho l} t\\
&=&\frac{8 \Delta p}{\rho l\pi^2}\sum_{n=1}^\infty \left( \frac{1}{n^2} – \frac{1}{(2n)^2} \right)\\
&=&\frac{8 \Delta p}{\rho l\pi^2} \frac{3}{4} \zeta(2)\\
&=&\frac{\Delta p}{\rho l} t
\end{eqnarray}
$$
となる。(無限級数の計算にはリーマンゼータ関数の特殊値\(\zeta(2)=\pi^2/6\)を用いた。)また、近似解は
$$
\bar{u}_{\rm app}(t)\sim C_{\rm cor} \frac{\Delta p}{\rho l} t
$$
となるため、比較すると\(C_{\rm cor}=1\)となる。

まとめると、過渡現象直近における現象を評価したい場合は補正係数なし(\(C_{\rm cor}=1\))で近似すべきであり、定常状態になるまでの現象まで考えたい場合は補正係数は\(C_{\rm cor}=0.834(1)\)を用いた近似を行うと良いことがわかった。

ポワズイユ流の非定常解


poiseuille
前回の定常解の導出において、解は初等的な積分で導出できた。しかし非定常解を求めるには一工夫が必要だ。
$$\frac{\partial u}{\partial t} = \frac{\Delta p}{\rho l} + \frac{\mu}{\rho} \frac{\partial^2 u}{\partial y^2}$$
の右辺第一項がない斉次方程式は、拡散方程式と呼ばれており、Fourier変換を行うか、あるいは変数分離解を仮定することなどにより求められる。

斉次方程式の解の導出

では、まずは斉次方程式

$$\frac{\partial \tilde{u}}{\partial t} – \frac{\mu}{\rho}\frac{\partial^2 \tilde{u}}{\partial y^2} = 0 $$

の解を求める。\( \tilde{u}(y,t) := \tilde{Y}(y) \tilde{T}(t) \)とし、変数分離解を求める。

$$\frac{1}{\tilde{T}}\frac{d \tilde{u}}{d t} = \frac{\mu}{\rho}\frac{1}{\tilde{Y}}\frac{d^2 \tilde{Y}(y)}{dy^2} = -C_0$$

とする(境界条件を満たす解が存在するよう、符号を選んだ。)と、

$$\begin{eqnarray}
\begin{cases}
\frac{\tilde{T}(t)}{dt}=-C_0\tilde{T}(t) \Longrightarrow \tilde{T}(t) \propto \exp (-C_0 t)\\
\frac{d^2 \tilde{Y}(y)}{dy^2} = -C_0 \frac{\rho}{\mu}\tilde{Y}(y) \Longrightarrow \tilde{Y}(y) \propto \sin \left(\sqrt{C_0\frac{\rho}{\mu}}y\right)
\end{cases}
\end{eqnarray}$$

と表される。

非斉次項による境界条件の変化

通常、ここで境界条件を課すことで解を満足する\(C_0\)を求めることができる。しかし今回の場合は斉次方程式の解を求める際に右辺第一項があるせいで境界条件を適用できず、\(C_0\)の条件を求められない。
そこで、右辺第一項は\(0 < y < d\)においてのみその値を持ち、\(y=0,\ d\)においてはゼロとなる関数へと変えてしまえば\(y=0, d\)で\(\tilde{u}(y,t)=0\)という境界条件を、斉次方程式に対して与えられる。

すると、\(\sqrt{C_0\frac{\rho}{\mu}}=\frac{n\pi}{d}\)となる。ところで、今回は境界条件より、解は\(y=0, d\)においてゼロとなり、かつ、与式は\(y=1/2\)において対称となっている。つまり、斉次方程式を解を導出するという議論は不要であり、解はFourier級数展開を用いて\( u(y,t)=\sum_{n=1}^\infty T_n(t)\sin \left(\frac{n\pi}{d}y\right) \)の形で表されることは明らかである。

境界条件を満たせる関数は右辺第一項は展開係数\(f_n\)を用いて
$$ \frac{\Delta p}{\rho l} = \sum_{n=1}^\infty f_n \sin\left(\frac{n\pi}{d}y\right)$$
とFourier級数展開で表されることがわかる。展開係数は\( \sin\left(\frac{m\pi}{d}y\right) \)をかけて\(0\rightarrow d\)の範囲で積分することで求められ、

$$f_n=\frac{2}{m\pi}\frac{\Delta p}{\rho l}\{1-(-1)^n \}$$

となる。最後に、これらを非斉次方程式に代入して、

$$
\sum_{n=1}^\infty \frac{\partial T_n(t)}{\partial t} \sin\left(\frac{n\pi}{d}y\right) = \sum_{n=1}^\infty f_n \sin\left(\frac{n\pi}{d}y\right) – \sum_{n=1}^\infty \frac{\mu}{\rho}\frac{n^2\pi^2}{d^2} T_n(t) \sin\left(\frac{n\pi}{d}y\right)
$$
直交性より、
$$
\frac{\partial T_n(t)}{\partial t} = f_n – \frac{\mu}{\rho}\frac{n^2\pi^2}{d^2} T_n(t)
$$
右辺全体を\(\tau_n\)と変数変換すれば
$$
\frac{\partial \tau_n(t)}{\partial t}=-\frac{\mu n^2\pi^2}{\rho d^2}\tau_n(t)
$$
となり、
$$
\tau_n(t)=C_{1n} \exp\left(-\frac{\mu n^2\pi^2}{\rho d^2} t \right)
$$
で表される。\(t=0\)において\(T_n=0\)とすると、積分定数は\(C_{1n}= f_n\)と決定され、\(T_n(t)\)は
$$
T_n(t)=f_n\frac{\rho d^2}{\mu n^2\pi^2}\left( 1- \exp\left( -\frac{\mu n^2\pi^2}{\rho d^2} t \right) \right)
$$
となる。まとめると、一般解は
$$
u(y,t)=\sum_{n=1}^\infty f_n\frac{\rho d^2}{\mu n^2\pi^2}\left( 1- \exp\left( -\frac{\mu n^2\pi^2}{\rho d^2} t \right) \right)\sin \left(\frac{n\pi}{d}y\right)
$$
で表される。最後に、\(f_n\)は、\(n\)が奇数の場合のみ値を持つ。そのため、整数\(m\)を用いて\(n=2m-1\)と置き直し、\(m\)を\(n\)で書き直すと
$$
u(y,t)=\sum_{n=1}^\infty \frac{4 d^2\Delta p}{\mu l(2n-1)^3\pi^3} \left( 1- \exp\left( -\frac{\mu (2n-1)^2\pi^2}{\rho d^2} t \right) \right)\sin \left(\frac{(2n-1)\pi}{d}y\right)
$$
となる。

検証

定常解との比較

今回の解の正しさを検証する。まず、\(t\to\infty\)としたときに定常解と一致するかを調べる。前回の導出の結果を、Fourier級数展開して比較を行う。
$$u(y)=-\frac{\Delta p}{2\mu l}y(y-d)=\sum_{n=1}^\infty g_n\sin \left(\frac{(2n-1)\pi}{d}y\right)$$
展開係数\(g_n\)は、
$$
\begin{eqnarray}
g_n&=&-\frac{2}{d}\int_0^d \frac{\Delta p}{2\mu l}y(y-d)\sin \left(\frac{(2n-1)\pi}{d}y\right)dy\\
&=&\frac{4 d^2 \Delta p }{\mu l (2 n-1)^3\pi^3}
\end{eqnarray}
$$
これは、今回の導出結果
$$
u(y,t\to\infty)=\sum_{n=1}^\infty \frac{4 d^2\Delta p}{\mu l(2n-1)^3\pi^3} \sin \left(\frac{(2n-1)\pi}{d}y\right)
$$
と一致していることが確認できた。

時間発展

この解はどのような時間発展をするのかを観察するため、平均速度\(\bar{u}(t):=\frac{1}{d}\int_0^d u(y,t) dy\)を求めると、
$$
\bar{u}(t)=\sum_{n=1}^\infty \frac{8 d^2\Delta p}{\mu l(2n-1)^4\pi^4} \left( 1- \exp\left( -\frac{\mu (2n-1)^2\pi^2}{\rho d^2} t \right) \right)
$$
となる。

ポワズイユ流の定常解


poiseuille

このように、長さ\(l\)、幅\(d\)の二次元平板間に圧力差\(\Delta p := p_{\rm in} – p_{\rm out}\)があるポワズイユ流の定常解を求めてみよう。
まず、ナビエ・ストークス方程式は非圧縮であるとき、速度\(\boldsymbol{v}\)は密度\(\rho\)、圧力\(p\)、粘性\(\mu\)、体積力\(\boldsymbol{K}\)を用いて
$$\frac{D\boldsymbol{v}}{Dt} =-\frac{1}{\rho} \nabla p + \frac{\mu}{\rho}\nabla^2 \boldsymbol{v} + \boldsymbol{K}$$
で表される。今回、\(\nabla p = -\frac{\Delta p}{l}\)として\(x\)成分を書き下すと、
$$\frac{\partial u}{\partial t} = \frac{\Delta p}{\rho l} + \frac{\mu}{\rho} \frac{\partial^2 u}{\partial y^2}$$
となる。ただし速度の成分は\(\boldsymbol{v}(u,v)\)とした。

定常解を求めるためには、まず左辺を\(0\)として考えればよく、
$$\frac{d^2 u}{d y^2} =-\frac{\Delta p}{\mu l}$$
と簡単な全微分で表せる。よって、積分して境界条件を考慮すれば、
$$u(y)=-\frac{\Delta p}{2\mu l}y(y-d)$$
となることがわかり、二次関数的な速度分布であるといえる。