(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)\)を用いた近似を行うと良いことがわかった。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

日本語が含まれない投稿は無視されますのでご注意ください。(スパム対策)