多重入れ子メッシュを使った解析

By M.Sato

August 1, 2024

概要

流体の内部に回転する物体がある場合を、Ansys Fluentに適用することを考える。

以前紹介したブログ記事においては、メッシュが入れ子関係になっており、親メッシュが回転しつつ、同時に子メッシュも回転している場合を想定し解析の手順を示した。

これは入れ子構造が一個しか存在しない最も単純な場合である。

今回は、入れ子構造を一般化し多重入れ子になったメッシュがそれぞれ独立に回転移動するような場合を紹介する。

例えば、太陽の周りを地球が回転移動し、地球の周囲を月が回転移動し、月の周囲を周回人工衛星が周回しているような場面を想像してもらえればよい。

ここでは三重入れ子の場合を紹介しているが一般の多重入れ子に対しても容易に拡張可能である。

メッシュの運動の記述

ここでは、領域が"out",“mid”,“in"の三個に分かれており、それぞれが独立に回転している状況を考える。

回転モデル

位置ベクトル

領域"out”

半径$\rho$、時刻$t$における座標点$\mathbf{r}_{out}(\rho,t)$は次式で表される。

ここで、$\mathbf{R}_{out}(t)$は領域"out"の回転中心位置ベクトルであり、移動しないものとする。

$\omega_0$は、回転角速度、$\theta_0$は初期回転角である。角度の基準軸は$+x$軸方向とする。

$$ \begin{align} \mathbf r_{out}(\rho,t) &= \mathbf R_{out}(t)+ \begin{pmatrix} \rho \cdot cos(\theta) \\ \rho \cdot sin(\theta) \end{pmatrix} \\ &= \begin{pmatrix} x_0 \\ y_0 \end{pmatrix} + \begin{pmatrix} \rho \cdot cos(\omega_0 t + \theta_0) \\ \rho \cdot sin(\omega_0 t + \theta_0) \end{pmatrix} \end{align} $$

領域"mid"

次に領域"mid"の座標点を表す式を導く。

半径$\rho$、時刻$t$における座標点$\mathbf{r}_{mid}(\rho,t)$は次式で表される。

領域"mid"は、回転中心$\mathbf{R}_{mid}(t)$を中心にして角速度$\omega_1$で回転している。

そして$\mathbf R_{mid}(t)$は、$\mathbf{R}_{out}(t)$の周りを角速度$\omega_0$で回転している。

$\mathbf R_{mid}(t)$と$\mathbf R_{out}(t)$の差分ベクトルは$\mathbf a_1$で表される。

静止座標系から見た、領域"mid"の回転角速度は$\omega_0+\omega_1$となる。 $\theta_1$は領域"mid"の初期回転角である。

$$ \begin{align} \mathbf v_{mid}(\rho,t) = \frac{d}{dt}\left(\mathbf r_{mid}(\rho,t)\right) &= \begin{pmatrix} -a_1 \omega_0 \cdot sin(\omega_0 t + \theta_0) \\ a_1 \omega_0 \cdot cos(\omega_0 t + \theta_0) \end{pmatrix} \\ & + \begin{pmatrix} -\rho \omega_1 \cdot sin(\omega_1 t + \theta_1) \\ \rho \omega_1 \cdot cos(\omega_1 t + \theta_1) \end{pmatrix} \\ \end{align} $$

領域"in"

$$ \begin{align} \mathbf v_{in}(\rho,t) = \frac{d}{dt}\left(\mathbf r_{in}(\rho,t)\right) &= \begin{pmatrix} -a_1 \omega_0 \cdot sin(\omega_0 t + \theta_0) \\ a_1 \omega_0 \cdot cos(\omega_0 t + \theta_0) \end{pmatrix} \\ & + \begin{pmatrix} -a_2 \omega_1 \cdot sin(\omega_1 t + \theta_1) \\ a_2 \omega_1 \cdot cos(\omega_1 t + \theta_1) \end{pmatrix} \\ & + \begin{pmatrix} -\rho \omega_2 \cdot sin(\omega_2 t + \theta_2) \\ \rho \omega_2 \cdot cos(\omega_2 t + \theta_2) \end{pmatrix} \\ \end{align} $$

領域"out", “mid”, “in” の各速度ベクトルの式を見ると最後の項

$$ (-\rho\omega\cdot sin(\omega t+\theta), \;\rho\omega\cdot cos(\omega t+\theta)) $$

は回転円板自体の回転速度を表しており、これ以外の項は多重入れ子回転円板の中心点の速度ベクトルを表している。

fluentにおける解析設定(単純な入れ子)

まず不連続インターフェイスの指定を行う。

次に単純な親子関係（１重入れ子）の場合、fluentのセル設定において以下のような相対指定を行うことが可能である。

ここではセルゾーン基準に親メッシュゾーンを指定している。

回転モデル

残念ながら多重入れ子を使う場合、fluentが未対応のためここで示したGUI指定は使えない。

代わりに以下に示すようにUDF（ユーザ定義関数）を使って指定を行う。

fluentにおける解析設定(多重入れ子)

インターフェイスの指定

不連続インターフェイスの指定は従来と同様の手順で行う。

インターフェース

セルゾーンの指定

ゾーン運動関数に各ゾーンに対応したUDF関数名を指定する。

今回の解析領域は三個のセルゾーンから構成されていることから、三種類のUDF関数を指定する。

fluid-out

fluid-mid

fluid-in

UDFの記述

セルゾーンの移動を定義するために"DEFINE_ZONE_MOTION"という関数マクロを使用する。

使用法

DEFINE_ZONE_MOTION (name, omega, axis, origin, velocity, time, dtime)

引数の型	説明
symbol name	UDFの名前．
real *omega	回転速度の大きさを指すポインタ．デフォルトは$0$．
real axis[3]	二次元軸対象の回転軸方向のベクトル．デフォルトは$(0\;0\;1)$と$(1\;0\;0)$．
real origin[3]	回転軸原点ベクトル．デフォルトは$(0\;0\;0)$．
real velocity[3]	平行移動速度ベクトル．デフォルトは$(0\;0\;0)$．
real current_time	現在時間．
real dtime	現在の時間ステップ．

以下に示す解析例での記述例を示す。

originおよびvelocityは各ゾーンの回転中心の挙動を記述していることに注意。

ARM1は “out"と"mid"の中心点間の距離、ARM2は “mid"と"in"の中心点間の距離、を表している。

OMG0,OMG1,OMG2はそれぞれ"out”、“mid”、“in"の相対回転角速度である(すなわち親ゾーンから見たときの子ゾーンの回転角速度)。

なお以下の例題では上記の位置ベクトル、速度ベクトルの説明とは異なり。角度の基準軸を$+y$軸方向とした記述になっているので注意。

#include "udf.h"

#define ARM1  10
#define ARM2  10
#define OMG0  (1.0)
#define OMG1  (2.0)
#define OMG2  (3.0)

/************************************************************/
DEFINE_ZONE_MOTION(out,omega,axis,origin,velocity,time,dtime)
/************************************************************/
{
   *omega = OMG0;

   axis[0] = 0.0;
   axis[1] = 0.0;
   axis[2] = 1.0;

   origin[0] = 0.0;
   origin[1] = 0.0;
   origin[2] = 0.0;

   velocity[0] = 0;
   velocity[1] = 0;
   velocity[2] = 0;

   return;
}


/************************************************************/
DEFINE_ZONE_MOTION(mid,omega,axis,origin,velocity,time,dtime)
/************************************************************/
{
   real theta0;

   *omega = OMG0 + OMG1;
   theta0 = OMG0 * time;

   axis[0] = 0.0;
   axis[1] = 0.0;
   axis[2] = 1.0;

   origin[0] = -ARM1 * sin(theta0);
   origin[1] = +ARM1 * cos(theta0);
   origin[2] =  0.0;

   velocity[0] = -ARM1 * OMG0 * cos(theta0);
   velocity[1] = -ARM1 * OMG0 * sin(theta0);
   velocity[2] = 0;

   return;
}


/************************************************************/
DEFINE_ZONE_MOTION(in,omega,axis,origin,velocity,time,dtime)
/************************************************************/
{
   real theta0,theta1;

   *omega = OMG0 + OMG1 + OMG2;
   theta0 = OMG0 * time;
   theta1 = (OMG0+OMG1) * time;

   axis[0] = 0.0;
   axis[1] = 0.0;
   axis[2] = 1.0;

   origin[0] = -ARM1 * sin(theta0)
               -ARM2 * sin(theta1);
   origin[1] = +ARM1 * cos(theta0)
               +ARM2 * cos(theta1);
   origin[2] = 0.0;

   velocity[0] = -ARM1 * OMG0 * cos(theta0)
                 -ARM2 * (OMG0+OMG1) * cos(theta1);
   velocity[1] = -ARM1 * OMG0 * sin(theta0)
                 -ARM2 * (OMG0+OMG1) * sin(theta1);
   velocity[2] = 0;

   return;
}