月別アーカイブ: 2022年9月

渦度法


\(\def\bm{\boldsymbol}\)\(\def\di{\displaystyle}\)\(\def\ve{\varepsilon_0}\)\(\newcommand{\pdr}[2]{\dfrac{\partial #1}{\partial #2}}\)\(\newcommand{\ppdr}[2]{\dfrac{\partial^2 #1}{\partial #2}}\)
◆ はじめに
非圧縮性の流れ方程式は、連続の式とナビエストークス方程式(無次元化した)を連立して解くことになる。
\(\left\{\begin{array}{l}\nabla\cdot\bm{V}=0\hspace{53mm}\cdots\,(1)\\
\pdr{\bm{V}}{t}+(\bm{V}\cdot\nabla)\bm{V}=-\nabla p+\dfrac{1}{\mathrm{Re}}\Delta\bm{V}\qquad\cdots\,(2)\end{array}\right.\)

この方程式は非線形で、かつ速度\(\,\bm{V}\,\)については時間発展形になっているが、圧力\(\,p\,\)についてはそうなっていない。そこで、速度\(\,\bm{V}\,\)を時間発展的に求める場合、各時間ステップでの連続の式(式(1))を満足するように圧力\(\,p\,\)を求める必要がある。
非圧縮性ナビエストークス方程式を数値的に解く方法はいくつかあるが、ここでは圧力を消去する渦度法をここで説明する。

◆ 導出
式(1)(2)を2次元で表すと
\(\left\{\begin{array}{l}\pdr{u}{x}+\pdr{v}{y}=0\hspace{40mm}\cdots\,(3)\\
\pdr{u}{t}+u\pdr{u}{x}+v\pdr{u}{y}=-\pdr{p}{x}+\dfrac{1}{\mathrm{Re}}\left(\ppdr{u}{x^2}+\ppdr{u}{y^2}\right)\qquad\cdots\,(4)\\
\pdr{v}{t}+u\pdr{v}{x}+v\pdr{v}{y}=-\pdr{p}{y}+\dfrac{1}{\mathrm{Re}}\left(\ppdr{v}{x^2}+\ppdr{v}{y^2}\right)\qquad\,\cdots\,(5)\end{array}\right.\)

まず、式(4)を\(\,y\,\)で偏微分する。左辺は
\(\quad\pdr{}{t}\left(\pdr{u}{y}\right)+\pdr{u^2}{x\,\partial y}+u\,\ppdr{u}{x\,\partial y}+\pdr{v}{y}\pdr{u}{y}+v\,\ppdr{u}{y^2}\qquad\)となる。変形すると
\(\quad=\pdr{}{t}\left(\pdr{u}{y}\right)+u\,\ppdr{u}{x\,\partial y}+\pdr{u}{y}\left(\pdr{u}{x}+\pdr{v}{y}\right)+v\,\ppdr{u}{y^2}\qquad\)式(3)を代入すると
\(\quad=\pdr{}{t}\left(\pdr{u}{y}\right)+u\,\ppdr{u}{x\,\partial y}+v\,\ppdr{u}{y^2}\qquad\cdots\,(6)\quad\)となる。
右辺は
\(\quad -\ppdr{p}{x\,\partial y}+\dfrac{1}{\mathrm{Re}}\,\pdr{}{y}\left(\ppdr{u}{x^2}+\ppdr{u}{y^2}\right)\qquad\cdots\,(7)\quad\)となる。

続けて、式(5)を\(\,x\,\)で偏微分する。左辺は式(4)と同様に変形すると
\(\quad\pdr{}{t}\left(\pdr{v}{x}\right)+u\,\ppdr{v}{x^2}+v\,\ppdr{v}{x\,\partial y}\qquad\cdots\,(8)\quad\)となる。
右辺は
\(\quad -\ppdr{p}{x\,\partial y}+\dfrac{1}{\mathrm{Re}}\,\pdr{}{x}\left(\ppdr{v}{x^2}+\ppdr{v}{y^2}\right)\qquad\cdots\,(9)\quad\)となる。

ここで、式(8)から式(6)を引いて偏微分の順序を入れ替えると
\(\quad\pdr{}{t}\left(\pdr{v}{x}-\pdr{u}{y}\right)+u\,\pdr{}{x}\left(\pdr{v}{x}-\pdr{u}{y}\right)+v\,\pdr{}{y}\left(\pdr{v}{x}-\pdr{u}{y}\right)\qquad\cdots\,(10)\quad\)となる。

式(9)から式(7)を引くと
\(\quad\dfrac{1}{\mathrm{Re}}\left(\pdr{}{x}\ppdr{v}{x^2}+\pdr{}{x}\ppdr{v}{y^2}-\pdr{}{y}\ppdr{u}{x^2}-\pdr{}{y}\ppdr{u}{y^2}\right)\)
\(\quad=\dfrac{1}{\mathrm{Re}}\left(\ppdr{}{x^2}\pdr{v}{x}+\ppdr{}{y^2}\pdr{v}{x}-\ppdr{}{x^2}\pdr{u}{y}-\ppdr{}{y^2}\pdr{u}{y}\right)\)
\(\quad=\dfrac{1}{\mathrm{Re}}\left[\ppdr{}{x^2}\left(\pdr{v}{x}-\pdr{u}{y}\right)+\ppdr{}{y^2}\left(\pdr{v}{x}-\pdr{u}{y}\right)\right]\qquad\cdots\,(11)\quad\)となる。

ここで渦度\(\,\omega\,\)を流速ベクトル\(\,\bm{V}\,\)の回転とする。
\(\quad\omega=\nabla\times\bm{V}=\pdr{v}{x}-\pdr{u}{y}\hspace{40mm}\cdots\,(12)\)
この式(12)を式(10)と式(11)に代入して等号で結ぶと
\(\quad\pdr{\omega}{t}+u\,\pdr{\omega}{x}+v\,\pdr{\omega}{y}=\dfrac{1}{\mathrm{Re}}\left(\ppdr{\omega}{x^2}+\ppdr{\omega}{y^2}\right)\qquad\cdots\,(13)\)

流れ関数\(\,\phi(x,y)\,\)を考える。流速ベクトルを流れ関数の回転と定義する。
\(\quad\bm{V}=\nabla\times\,(0,0,\phi)=\left(\pdr{\phi}{y}\,,\,-\pdr{\phi}{x}\,,\,0\right)\qquad\cdots\,(14)\)
2次元の場合、\(\bm{V}=(\,u\,,\,v\,,\,0)\,\)の渦度の\(\,z\,\)成分\(\,\omega\,\)は
\(\quad\omega=\pdr{v}{x}-\pdr{u}{y}=-\left(\ppdr{\phi}{x^2}+\ppdr{\phi}{y^2}\right)\qquad\cdots\,(15)\)

よって、式(1)(2)は流れ関数\(\,\phi\,\)と渦度\(\,\omega\,\)の連立方程式として、以下の形になる。
\(\left\{\begin{array}{l}\ppdr{\phi}{x^2}+\ppdr{\phi}{y^2}=-\omega\hspace{60mm}\cdots\,(16)\\
\pdr{\omega}{t}+\pdr{\phi}{y}\,\pdr{\omega}{x}-\pdr{\phi}{x}\,\pdr{\omega}{y}=\dfrac{1}{\mathrm{Re}}\left(\ppdr{\omega}{x^2}+\ppdr{\omega}{y^2}\right)\qquad\cdots\,(17)\end{array}\right.\)

◆ 解き方
式(16)、(17)を差分近似する。差分については「差分について」を参照のこと。

\(\quad a_2\phi_{i-1,j}+b_2\phi_{i,j}+c_2\phi_{i+1,j}+a_4\phi_{i,j-1}+b_4\phi{i,j}+c_4\phi_{i,j+1}=-\omega_{i,j}\qquad\cdots\,(18)\)
\(\quad\dfrac{\omega_{i,j}^{n+1}-\omega_{i,j}}{\Delta t}+(a_3\phi_{i,j-1}+b_3\phi_{i,j}+c_3\phi_{i,j+1})(a_1\omega_{i-1,j}+b_1\omega_{i,j}+c_1\omega_{i+1,j})\)
\(\hspace{28mm}-(a_1\phi_{i-1,j}+b_1\phi{i,j}+c_1\phi_{i+1,j})(a_3\omega_{i,j-1}+b_3\omega_{i,j}+c_3\omega_{i,j+1})\)
\(\hspace{15mm}=\dfrac{1}{\mathrm{Re}}(a_2\omega_{i-1,j}+b_2\omega_{i,j}+c_2\omega_{i+1,j}+a_4\omega_{i,j-1}+b_4\omega_{i,j}+c_4\omega_{i,j+1})\qquad\cdots\,(19)\)

まず、式(18)をSOR法(参照)で解くために、式(18)を\(\,\phi_{i,j}\,\)について解いて

\(\quad\phi_{i,j}^*=-\dfrac{1}{b_2+b_4}(a_2\phi_{i-1,j}^{(\nu+1)}+c_2\phi_{i+1,j}^{(\nu)}+a_4\phi_{i,j-1}^{(\nu+1)}+c_4\phi_{i,j+1}^{(\nu)}+\omega_{i,j})\qquad\cdots\,(20)\)

\(\quad\phi_{i,j}^{(\nu+1)}=(1-\alpha)\phi_{i,j}^{(\nu)}+\alpha\phi_{i,j}^*\qquad\cdots\,\,(21)\)

として、反復解法(SOR法)で求める。ただし、\(\,\nu\,\)は反復回数、\(\,\alpha\,\)は加速係数である。


参考文献

  • 流体解析の基礎:河村 哲也(2014)、朝倉書店
  • 流体力学第2版:杉山 弘、遠藤 剛、新井 隆景(2014)、森北出版

流れ関係特性値


\(\def\bm{\boldsymbol}\)\(\def\di{\displaystyle}\)\(\def\ve{\varepsilon_0}\)\(\newcommand{\pdr}[2]{\dfrac{\partial #1}{\partial #2}}\)\(\newcommand{\ppdr}[2]{\dfrac{\partial^2 #1}{\partial #2}}\)
◆ レイノルズ数 ( Re )
レイノルズ数(Reynolds number)は流体力学において慣性力と粘性力との比で定義される無次元量である。流れの中でのこれら2つの力の相対的な重要性を定量している。

\(\quad\mathrm{Re}=\dfrac{\rho\,U\,L}{\mu}\)

  • \(\rho\):流体の密度 kg/m\(^3\quad\) (例:水:\(1000\) kg/m\(^3\) 空気:\(1.2\) kg/m\(^3\))
  • \(U\):流体の速度 m/s
  • \(L\):代表長さ m
  • \(\mu\):流体の粘度 Pa\(\cdot\)s , kg/m\(\cdot\)s\(\quad\) (例:水:\(1.01\) mPa\(\cdot\)s 空気:\(0.0181\) mPa\(\cdot\)s at \(20\)℃)
  • \(\star\)
    【例】内径 \(2.5\) cm の配管に水が流れた場合、流速 \(U=0.04\) m/s の時は
    \(\quad \mathrm{Re}=\dfrac{1000\times 0.04\times 0.025}{1.01\times 10^{-3}}=1000\quad\)となる。
    \(\qquad\)空気の場合、Re=\(1000\) とするには、流速は \(U=0.6\) m/s 必要となる。
    \(\quad \mathrm{Re}=\dfrac{1.2\times 0.6\times 0.025}{0.018\times 10^{-3}}=1000\quad\)となる。

    ◆ プラントル数 ( Pr )
    プラントル数( Prandtl number)は熱伝導に関する無次元の物性値であり、熱エネルギー拡散の度合いと運動エネルギー拡散の度合いの比です。プラントル数が大きいほど熱エネルギーの拡散に対して運動エネルギーの拡散の度合いが強いことを意味します。プラントル数は流体に固有の物性値です。流動現象において見られる壁面近傍の境界層を例にすると、プラントル数が \(1.0\) の場合、温度境界層と速度境界層の厚みは等しくなります。プラントル数が \(1.0\) より大きくなると、温度境界層に対して速度境界層の方が厚くなり、プラントル数が \(1.0\) より小さい場合はその逆となります。

    \(\quad\mathrm{Pr}=\dfrac{\nu}{\alpha}=\dfrac{\mu\,C_P}{k}\)

  • \(\nu=\mu/\rho\):流体の動粘度 m\(^2\)/s
  • \(\alpha\):熱拡散率(温度伝導率) m\(^2\)/s
  • \(\mu\):流体の粘度 Pa\(\cdot\)s , kg/m\(\cdot\)s\(\quad\) (例:水:\(1.01\) mPa\(\cdot\)s 空気:\(0.0181\) mPa\(\cdot\)s at \(20\)℃)
  • \(C_P\):流体の比熱(定圧比熱) J/kg\(\cdot\)K\(\quad\)(例:水:\(4182\) J/kg\(\cdot\)K 空気:\(1006\) J/kg\(\cdot\)K at \(20\)℃)
  • \(k\):熱伝導率 W/m\(\cdot\)K , J/s\(\cdot\)m\(\cdot\)K\(\quad\)(例:水:\(0.602\) W/m\(\cdot\)K 空気:\(0.0257\) W/m\(\cdot\)K at \(20\)℃)
  • \(\star\)
    水(\(300\)K \(1\)atm):\(5.466\) 空気(\(300\)K \(1\)atm):\(0.717\) 水銀(\(300\)K \(1\)atm):\(0.025\) 潤滑油(\(300\)K \(1\)atm):\(2040\)

    ◆ グランホフ数 ( Gr )
    グラスホフ数(Grashof Number)は、伝熱現象、物質移動現象に関して、流れ場における粘性力に対する浮力の相対的な影響を示す無次元量である。粘性力と浮力の大きさの比を表した無次元数です。自然対流において流れが層流から乱流に遷移するかどうかを特徴付けます。レイノルズ数が強制対流の流れを特徴付ける無次元数ならば、グラスホフ数は自然対流の流れを特徴付ける無次元数です。グラスホフ数が大きい、つまり浮力の影響が強い場合は、流体はより自由に流れようとするため流動は乱流場となります。一方、グラスホフ数が小さい場合は、流体の粘度による流れの抑制効果が高いため層流場となります。層流と乱流の境界となるグラスホフ数はアプリケーションによって異なりますが、概ね \(10^8\)~\(10^9\) のオーダーです。

    \(\quad\mathrm{Gr}=\dfrac{g\beta\Delta\theta L^3}{\nu^2}=\dfrac{g\rho^2\beta\Delta\theta L^3}{\mu^2}\)

  • \(g\):重力加速度 m/s\(^2\) \(9.807\) m/s\(^2\)
  • \(\beta\):体膨張係数 \(1\)/K (例:水:\(0.2\times 10^{-3}/\)K 空気:1/K)
  • \(\Delta\theta\):温度差 K
  • \(\star\)
    空気で代表温度差:\(100\) K 代表長さ:\(1\) m では \(\mathrm{Gr}=1.466\times 10^{10}\) となる。

    差分について


    \(\def\bm{\boldsymbol}\)\(\def\di{\displaystyle}\)\(\def\ve{\varepsilon_0}\)\(\newcommand{\pdr}[2]{\dfrac{\partial #1}{\partial #2}}\)\(\newcommand{\ppdr}[2]{\dfrac{\partial^2 #1}{\partial #2}}\)
    ◆ はじめに
    流れ解析において、微分方程式の数値解を求める手法として、まず『差分法』を挙げることが出来る。 関数 \(f(x)\,\)の\(\,x=x_i\,\)における微分は

    \(\quad \dfrac{df}{dx}=\di\lim_{h\to 0}\dfrac{f(x_i+h)-f(x_i)}{h}\quad\cdots\,(1)\)

    で表される。ここで、極限をとらずに\(\,h\,\)が十分小さいとして、”\(\,\sim\,\)” の記号を使って、以下のように表す。

    \(\quad \dfrac{df}{dx}\,\sim\,\dfrac{f(x_i+h)-f(x_i)}{h}\quad\cdots\,(2)\)

    図で言うと、点 \(P\,\)の傾きを\(\,PA\,\)の傾きとするのが『前進差分』である。逆に\(\,BP\,\)とするのが『後退差分』(式(3))である。

    \(\quad \dfrac{df}{dx}\,\sim\,\dfrac{f(x_i)-f(x_i-h)}{h}\quad\cdots\,(3)\)

    点\(\,P\,\)での微分を\(\,BA\,\)の傾きとするのが『中心差分』である。図の差を\(\,h=h_1=h_2\,\)とする。

    \(\quad \dfrac{df}{dx}\,\sim\,\dfrac{f(x_i+h)-f(x_i-h)}{2h}\quad\cdots\,(4)\)

    ◆ 2次微分
    式(4)を\(\,x\,\)で微分すると

    \(\quad \left.\dfrac{d^2f}{dx^2}\right|_{x=x_i}=\dfrac{df^{\prime}}{dx}\,\sim\,\dfrac{f^{\prime}(x_i+h)-f^{\prime}(x_i-h)}{h}\quad\cdots\,(5)\)

    ここで、\(f^{\prime}(x_i+h)\,\)を\(\,f(x_i)\,\)の前進差分(式(1))として展開すると、

    \(\quad \left.\dfrac{df}{dx}\right|_{x=x_i}\,\sim\,\dfrac{f(x_i+h)-f(x_i)}{h}\quad\cdots\,(6)\)

    \(f^{\prime}(x_i-h)\,\)を\(\,f(x_i)\,\)の後退差分(式(2))として展開すると、

    \(\quad \left.\dfrac{df}{dx}\right|_{x=x_i}\,\sim\,\dfrac{f(x_i)-f(x_i-h)}{h}\quad\cdots\,(7)\)

    式(6),(7)を式(5)に代入すると

    \(\quad \left.\dfrac{d^2f}{dx^2}\right|_{x=x_i}\,\sim\,\dfrac{f(x_i+h)-2f(x_i)+f(x_i-h)}{h^2}\quad\cdots\,(8)\)

    ◆ 不均等格子での差分
     左図のような\(\,h_1\ne h_2\,\)の場合、不均等格子となる。三点の線形結合で微分を差分で表すことを考える。\(\,f(x_i)\,\)を\(\,f_i\,\)と表すこととする。

    \(\,a\,f_{i-1}+b\,f_i+c\,f_{i+1}\quad\cdots\,(10)\)

    テーラー展開は次のようになる。

    \(\quad g(e)=g(d)+g^{\prime}(d)(e-d)+\dfrac{g^{\prime\prime}(d)}{2!}(e-d)^2+\cdots\)

    \(\,a\,f_{i-1}\,\)をテーラー展開すると

    \(\quad af(x_i-h_1)=af(x_i)+af^{\prime}(x_i)(-h_1)+a\dfrac{f^{\prime\prime}(x_i)}{2!}(-h_1)^2+\cdots\quad (11)\)

    \(\,c\,f_{i+1}\,\)をテーラー展開すると

    \(\quad cf(x_i+-h_2)=cf(x_i)+cf^{\prime}(x_i)(h_2)+c\dfrac{f^{\prime\prime}(x_i)}{2!}(h_2)^2+\cdots\quad (12)\)

    式(11)と(12)を式(10)に代入すると

    \(\quad a\,f_{i-1}+b\,f_i+c\,f_{i+1}=(a+b+c)f_i+(ch_2-ah_1)f_i^{\prime}+(ah_1^2+ch_2^2)\dfrac{f_i^{\prime\prime}}{2}+\cdots\,(13)\)

    これより、1階微分\(\,f_i^{\prime}\,\)を\(\,f_{i-1}\,\,,\,\,f_i\,\,,\,\,f_{i+1}\,\)の線形結合で表すには

    \(\,a+b+c=0\ \ ,\quad ch_2-ah_1=1\ \ ,\quad 1/2(ah_1^2+ch_2^2)=0\ \)とおいて、\(a\ ,\ b\ ,\ c\ \)を求める。

    \(\,\,a=-\dfrac{h_2}{h_1(h_1+h_2)}\,\,,\quad b=\dfrac{h_2-h_1}{h_1h_2}\,\,,\quad c=\dfrac{h_1}{h_2(h_1+h_2)}\quad\)となる。具体的な近似式は

    \(\quad\left.\dfrac{df}{dx}\right|_{x=x_i}=\dfrac{-h_2^2\,f_{i-1}+(h_2-h_1)(h_2+h_1)\,f_i+h_1^2\,f_{i+1}}{h_1h_2(h_1+h_2)}\quad\cdots\,(14)\)

    同様に、2階微分\(\,f_i^{\prime\prime}\,\)を\(\,f_{i-1}\,\,,\,\,f_i\,\,,\,\,f_{i+1}\,\)の線形結合で表すには

    \(\,a+b+c=0\ \ ,\quad ch_2-ah_1=0\ \ ,\quad 1/2(ah_1^2+ch_2^2)=1\ \)とおいて、\(a\ ,\ b\ ,\ c\ \)を求める。

    \(\,\,a=\dfrac{2}{h_1(h_1+h_2)}\,\,,\quad b=-\dfrac{2}{h_1h_2}\,\,,\quad c=\dfrac{2}{h_2(h_1+h_2)}\quad\)となる。近似式は以下のようになる。

    \(\quad\left.\dfrac{d^2f}{dx^2}\right|_{x=x_i}=\dfrac{2\,h_2\,f_{i-1}-2\,(h_2+h_1)\,f_i+2\,h_1\,f_{i+1}}{h_1h_2(h_1+h_2)}\quad\cdots\,(15)\)

    流れ解析の練習-2 室内換気 不均等格子


    \(\def\bm{\boldsymbol}\)\(\def\di{\displaystyle}\)\(\def\ve{\varepsilon_0}\)\(\newcommand{\pdr}[2]{\dfrac{\partial #1}{\partial #2}}\)\(\newcommand{\ppdr}[2]{\dfrac{\partial^2 #1}{\partial #2}}\)