渦輪の解析

By M.Sato

September 11, 2024

概要

ヘルムホルツの渦定理によれば、完全流体中に発生した渦には以下の特性がある。

渦糸に沿ってその強さ（渦度×渦糸断面積）は同じである
渦糸は流体中で消滅することはなく、境界部にまで伸びているか閉じたリング状になっているかである
渦のない運動は永遠に渦のないままで運動を続ける

リング状になった渦糸は渦輪とも呼ばれており、渦運動を利用して物質輸送を行うといった応用が考えられている。また空気砲といった科学おもちゃにも使われている

(参考：空気砲)

ここでは、Ansys Fluentを用いた数値計算により渦輪のシミュレーションを行ってみた。

解析の方針

渦輪は軸対称形状とみなすことができるため、軸対称モデルを使う
初期流れ場として、時刻ゼロにおいて瞬間的に渦輪が発生した場合を想定する
流体は空気を想定し、非圧縮、等温、層流を仮定する。

初期流れ場の生成

実際の渦輪発生装置は円形の流出口から流体を瞬間的に高速に噴出させるようになっているが、このような複雑な流れ場では渦輪単独での流動特性を把握するのが難しい。

従ってここでは、瞬間的に渦輪を発生させてその流れ場を初期条件として渦輪の挙動を計算した。

なお初期渦輪は非圧縮性流体の連続方程式を満足するように設定する（ソレノイダル場）。

ストークスの流れ関数

軸対称流れ場では流速をストークスの流れ関数$\Psi$を使って表すことができる。

簡単のため軸方向座標、流速を $x,u$、半径方向座標、流速を$y,v$で表すと流速は次式で表される。

$$ \begin{align} u &= \frac{1}{y}\frac{\partial \Psi}{\partial y} \\ v &= -\frac{1}{y}\frac{\partial \Psi}{\partial x} \end{align} $$

渦度$\omega$は次式で表される。

$$ \begin{align} \omega &= \frac{\partial v}{\partial x}-\frac{\partial u}{\partial y} \\ &= \frac{\partial}{\partial x}\left(-\frac{1}{y}\frac{\partial \Psi}{\partial x}\right)-\frac{\partial}{\partial y}\left(\frac{1}{y}\frac{\partial \Psi}{\partial y}\right) \end{align} $$

従って

$$ \begin{align} \frac{\partial}{\partial x}\left(\frac{1}{y}\frac{\partial \Psi}{\partial x}\right)+\frac{\partial}{\partial y}\left(\frac{1}{y}\frac{\partial \Psi}{\partial y}\right) &= -\omega \end{align} $$

これは、拡散係数を$1/y$にした二次元（軸対称ではない）ポアソン方程式になっている。

流れ関数の計算

$\Psi$を計算するために、まずは二次元モデルを作成した。

解析対象領域を以下に示す。 $Q_{\Psi}$と書かれた位置に渦管を配置する。

model-ini

メッシュ分割図を示す。全領域を１辺2mmの正方形で分割した。

mesh

$\Psi$に関するポアソン方程式をfluentを使って計算する。

$\omega$(渦度)は、渦輪の部分にだけ分布するようにした。

渦管は直径$0.03$[m]であり、渦度は渦管断面中心において最大値$1000$[1/s]となるように設定した。

また渦輪の位置は$x=y=0.07$[m]の位置とした。

渦度分布

メッシュ拡大

ユーザ定義関数の指定

$\Psi$方程式の拡散係数、ソース項をUDFを使い記述した。

また得られた $\Psi$場を空間微分することで流速を計算した。

なお$\Psi$変数はユーザ定義変数UDS0に割り当てている。

UDFのソースプログラム(fl-vort.c)を示す。

#include "udf.h"

#define  RADIMAX  0.015
#define  SRCMAX   1000.0
#define  XCENTER  0.07
#define  YCENTER  0.07


/*-----------------------------------------------------------------*/
DEFINE_DIFFUSIVITY(psidif,c,t,i)
/*-----------------------------------------------------------------*/
{
   real xyz[ND_ND];
   C_CENTROID(xyz, c, t);
   return 1.0/xyz[1];
}


/*-----------------------------------------------------------------*/
DEFINE_SOURCE(psisrc,c,t,dS,eqn)
/*-----------------------------------------------------------------*/
{
   real source,radi;
   real xyz[ND_ND];
   C_CENTROID(xyz, c, t);

   radi = sqrt((xyz[0]-XCENTER)*(xyz[0]-XCENTER)
              +(xyz[1]-YCENTER)*(xyz[1]-YCENTER));

   source = (radi>RADIMAX ? 0.0 : SRCMAX/RADIMAX*(RADIMAX-radi));
   dS[eqn] = 0.0;

   return source;
}


/*-----------------------------------------------------------------*/
DEFINE_ON_DEMAND(on_demand)
/*-----------------------------------------------------------------*/
{
   Domain *d = Get_Domain(1);
   Thread *t;
   cell_t c;
   real xyz[ND_ND];

   thread_loop_c (t,d)
   {
      if (NULL != THREAD_STORAGE(t,SV_UDS_I(0)) &&
          NULL != T_STORAGE_R_NV(t,SV_UDSI_G(0)))
      {
         begin_c_loop (c,t)
         {
            C_CENTROID(xyz, c, t);
            C_U(c,t) = C_UDSI_G(c,t,0)[1]/xyz[1];
            C_V(c,t) =-C_UDSI_G(c,t,0)[0]/xyz[1];
         }
         end_c_loop (c,t)
      }
   }

}