回転座標系を使った解析

By M.Sato

April 5, 2024

概要

回転する境界がある流体の流れを解析する場合、座標系を変換したほうが便利な場合がある。

例えば円筒形のタンク内を撹拌用のブレードが回転している場合、回転ブレードを静止させておき、周囲のタンクが回転していると考える方が簡単である。ただし、この手法はタンク内にバッフル板などの突起がなくタンクの運動を単純な回転境界として扱えることが必要である。

ここでは、回転座標系を使う場合の基礎方程式および解析例を紹介する。

ここで示す解析例は、単純な円筒形であるため、回転座標系を使わなくても解析可能であるが、検証のため、絶対（静止）座標系と回転座標系の場合で計算を行い妥当性を検証した。

回転座標系における運動

一般的な運動方程式の表現

ベクトル$\mathbf{a}$が回転軸周りを角速度$\Omega$で回転している場合を考える。$\mathbf{a}$の時間微分は次式で表される。

$$ \frac{d\mathbf{a}}{dt}=\mathbf{\Omega} \times \mathbf{a} \tag{1} $$

モデル

次に、任意のベクトル $\mathbf{f}$が静止座標系と回転座標系でどのように変換されるかを調べてみる。

$\mathbf{i},\mathbf{j},\mathbf{k}$を静止座標系の基底ベクトルとする。

$$ \mathbf{f}=f_1 \mathbf{i} + f_2 \mathbf{j} + f_3 \mathbf{k} \tag{2} $$ $$ \frac{d \mathbf{f}}{dt} = \frac{d f_1}{dt} \mathbf{i} + \frac{d f_2}{dt} \mathbf{j} + \frac{d f_3}{dt} \mathbf{k} + \frac{d \mathbf{i}}{dt} f_1 +\frac{d \mathbf{j}}{dt} f_2 + \frac{d \mathbf{k}}{dt} f_3 \notag $$ $$ = \frac{d f_1}{dt} \mathbf{i} + \frac{d f_2}{dt} \mathbf{j} + \frac{d f_3}{dt} \mathbf{k} + \Omega \times (f_1\mathbf{i}+f_2\mathbf{j}+f_3\mathbf{k}) \notag $$ $$ = \frac{d’ \mathbf{f}}{dt} +\Omega \times \mathbf{f} \tag{3} $$

なお、ここでは先述の $d\mathbf{i}/dt=\mathbf{\Omega} \times \mathbf{i}$などの関係式(例えば$\mathbf{a}=>\mathbf{i}$)を利用している。 $d’/dt$は回転座標系上での時間微分を表す。

$\Omega \times$を含めて微分作用素とみなすと上式は以下のように書くことができる。

$$ \frac{d \mathbf{f}}{dt} = \left( \frac{d’ }{dt} + \mathbf{\Omega} \times \right) \mathbf{f} \tag{4} $$

$$ \frac{d’ \mathbf{f}}{dt} = \left( \frac{d }{dt} - \mathbf{\Omega} \times \right) \mathbf{f} \tag{5} $$

$\mathbf{f}$を位置ベクトル$\mathbf{r}$とすると、次式が成り立つ。

$$ \frac{d \mathbf{r}}{dt} = \frac{d’ \mathbf{r}}{dt} + \mathbf{\Omega} \times \mathbf{r} \tag{6} $$ $$ \mathbf{v} = \mathbf{v}_R + \mathbf{\Omega} \times \mathbf{r} \tag{7} $$

ここで$\mathbf{v}$は静止座標系で表した速度、$\mathbf{v}_R$は回転座標系で表した速度である。

次に加速度を回転座標系で記述してみる。

$$ \frac{d^2\mathbf{r}}{dt^2} = \frac{d}{dt}\left( \frac{d \mathbf{r} }{dt} \right) = \frac{d’}{dt}\left( \frac{d’ \mathbf{r}}{dt} + \mathbf{\Omega} \times \mathbf{r}\right) +\mathbf{\Omega} \times \left( \frac{d’ \mathbf{r}}{dt} + \mathbf{\Omega} \times \mathbf{r}\right) \notag $$

$$ =\frac{d^{‘2} \mathbf{r}}{dt^2}+2\Omega\times \frac{d’\mathbf{r}}{dt} +\Omega\times\left(\Omega\times\mathbf{r}\right) +\frac{d\mathbf{\Omega}}{dt}\times \mathbf{r}　\notag $$

$$ =\frac{d^{‘2} \mathbf{r}}{dt^2}+2\Omega\times \mathbf{v}_R +\Omega\times\left(\Omega\times\mathbf{r}\right) +\dot{\mathbf{\Omega}} \times \mathbf{r} \tag{8} $$

これより、物体の運動方程式は次式で表される。

$$ m \frac{d^2\mathbf{r}}{dt^2} =m\frac{d^{‘2} \mathbf{r}}{dt^2} +2m\Omega\times \mathbf{v}_R +m\Omega\times\left(\Omega\times\mathbf{r}\right) +m\dot{\mathbf{\Omega}} \times \mathbf{r} = \mathbf{F} \tag{9} $$

$$ m\frac{d^{‘2} \mathbf{r}}{dt^2} = \mathbf{F} -2m\Omega\times \mathbf{v}_R -m\Omega\times\left(\Omega\times\mathbf{r}\right) -m\dot{\mathbf{\Omega}} \times \mathbf{r} \tag{10} $$

ここで $m$は物体の質量、$\mathbf{F}$は物体に作用する力である。
左辺は回転座標系で表した慣性力であり、右辺の第二項はコリオリ力、第三項は遠心力と呼ばれる項である。
第四項は回転角速度が時間変化する場合に現れる項である。

流体運動への適用

相対速度形式

これらの関係式を以下に示す連続方程式およびナビエ・ストークス式に適用することを考える。

$$ \nabla \mathbf{u} =0 \tag{11} $$ $$ \frac{D\mathbf{u}}{Dt} = -\frac{1}{\rho}\nabla p +\nu \nabla^2 \mathbf{u} +K \tag{12} $$

静止座標系の流速と回転座標系の流速には以下の関係式がある。
ここで $\mathbf{u}$は静止座標系での流速、$\mathbf{u}_R$は回転座標系での流速である。

$$ \mathbf{u} =\mathbf{u}_R + \Omega\times \mathbf{r} \tag{13} $$ $$ \mathbf{u}_R =\mathbf{u} - \Omega\times \mathbf{r} \tag{14} $$

連続方程式、運動方程式中に $\nabla \mathbf{u}$、および$\nabla^2 \mathbf{u}$が含まれるが、これらはそれぞれ以下のようになり$\Omega \times \mathbf{r}$に関する効果はゼロである。

$$ \nabla \mathbf{u} =\nabla \mathbf{u}_R+ \nabla \cdot( \Omega\times \mathbf{r} ) =\nabla \mathbf{u}_R \tag{15} $$

$$ \nabla^2 \mathbf{u} =\nabla^2 \mathbf{u}_R+ \nabla^2 ( \Omega\times \mathbf{r} ) =\nabla^2 \mathbf{u}_R \tag{16} $$

これらの関係式を使うと最終的な、連続方程式およびナビエ・ストークス式は以下のように表される。
ここで$\mathbf{K}$は外力項である。

$$ \nabla \mathbf{u}_R =0 \tag{17} $$

$$ \frac{D\mathbf{u}}{Dt} = \frac{D\mathbf{u}_R}{Dt}+2\Omega\times \mathbf{u}_R +\Omega\times\left(\Omega\times\mathbf{r}\right) +\dot{\mathbf{\Omega}} \times \mathbf{r} \notag $$

$$ =-\frac{1}{\rho}\nabla p +\nu \nabla^2 \mathbf{u} +\mathbf{K} \tag{18} $$

流速 $\mathbf{u}_R$に関する方程式として記述すると次式が得られる。

$$ \frac{D\mathbf{u}_R}{Dt} =\frac{\partial \mathbf{u}_R}{\partial t} +\nabla \left( \mathbf{u}_R \otimes \mathbf{u}_R \right) \notag $$ $$ = -\frac{1}{\rho}\nabla p +\nu \nabla^2 \mathbf{u}_R +\mathbf{K} -2\Omega\times \mathbf{u}_R -\Omega\times\left(\Omega\times\mathbf{r}\right) -\dot{\mathbf{\Omega}} \times \mathbf{r} \tag{19} $$

右辺の第四項はコリオリ力、第五項は遠心力と言われているものである。
この表現式は相対速度形式と呼ばれる。

絶対速度形式

次に絶対速度形式と呼ばれる定式化を導く。

$$ \frac{d^2\mathbf{r}}{dt^2} = \frac{d}{dt}\left( \frac{d \mathbf{r} }{dt} \right) =\frac{d}{dt}\left( \frac{d’ \mathbf{r} }{dt}+\Omega\times \mathbf{r} \right) \notag $$ $$ =\frac{d}{dt}\left(\frac{d’ \mathbf{r} }{dt} \right) +\frac{d \Omega}{dt}\times \mathbf{r} +\Omega \times \frac{d \mathbf{r} }{dt} \tag{20} $$

位置ベクトルの時間微分を流速で表し、オイラー座標で表現すると以下のようになる。

$$ \frac{D\mathbf{u}}{Dt} = \left[ \frac{\partial \mathbf{u}_R}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u}_R \right] +\frac{d\Omega}{dt}\times \mathbf{r} +\Omega \times \mathbf{u} \notag $$ $$ =\frac{\partial \mathbf{u}}{\partial t} -\frac{\partial }{\partial t}\left( \Omega \times \mathbf{r}\right) +\nabla \left( \mathbf{u} \otimes \mathbf{u}_R \right) +\frac{d\Omega}{dt}\times \mathbf{r} +\Omega \times \mathbf{u} \notag $$ $$ =\frac{\partial \mathbf{u}}{\partial t} +\nabla \left( \mathbf{u} \otimes \mathbf{u}_R \right) +\Omega \times \mathbf{u} \notag $$ $$ =-\frac{1}{\rho}\nabla p +\nu \nabla^2 \mathbf{u} +K \tag{21} $$

最終的に以下の絶対速度形式が得られる。

$$ \frac{\partial \mathbf{u}}{\partial t} +\nabla \left( \mathbf{u} \otimes \mathbf{u}_R \right) +\Omega \times \mathbf{u} =-\frac{1}{\rho}\nabla p +\nu \nabla^2 \mathbf{u} +K \tag{22} $$

解析例

ここでは fluentの検証例題#8 (Ansysユーザのみアクセス可) を対象として計算を行った。

元の例題データは軸対称モデルを使っているが、ここでは３次元モデル(1/4円筒)を使って解析を行ってみた。なおcase-srfでの流速は相対速度形式を使った。

解析モデル

解析モデル図を示す。

モデル

物性値

密度: $\rho=1.0[kg/m^3]$
粘度: $\mu=5.56\times 10^{-4}[Pa.s]$

メッシュ

1/4円筒の領域を解析領域としてメッシュを作成した。円筒切断面には周期境界条件を設定した。

メッシュ

境界条件

回転の効果を与えるために壁面に速度を指定する。
解析ケースにより異なる指定を行った。

検査

解析ケース

ここでは以下の３ケースで計算を行った。
SRFはSingle-Reference-Frameの意味である。

case	内容
nosrf	静止座標系での計算
srf	回転座標系での計算（SRF法使用）
srfudf	静止座標系での計算。ただしUDFで見かけの力を指定し、SRF法と同等の境界速度を指定

case-srfudfはfluentに組み込みのSRFオプションを使わないでUDFで体積力を指定して同等の計算ができるかを調べるものである。
計算に使用したUDFソースプログラムを以下に示す。

fl-srf.c:

#include "udf.h"

#define  OMEGA   1.0
#define  OMEGA2  (OMEGA*OMEGA)

/* origin of local coord */
#define  OX  0
#define  OY  0
#define  OZ  0

/* normalized omega axis */
#define  AX  0
#define  AY  0
#define  AZ  1

/*--------------------------------------------*/
DEFINE_SOURCE(xmom,c,t,dS,eqn)
/*--------------------------------------------*/
{
   real omgx,omgy,omgz;
   real relx,rely,relz;
   real velx,vely,velz;
   real colx,coly,colz;
   real cenx,ceny,cenz;
   real radx,rady,radz;
   real x[ND_ND],source,rxk;

   C_CENTROID(x,c,t);

   relx = x[0] - OX;
   rely = x[1] - OY;
   relz = x[2] - OZ;

   velx = C_U(c,t);
   vely = C_V(c,t);
   velz = C_W(c,t);

   omgx = OMEGA*AX;
   omgy = OMEGA*AY;
   omgz = OMEGA*AZ;

/* coriolis force */
   colx = 2.0*(omgy*velz - omgz*vely);
   coly = 2.0*(omgz*velx - omgx*velz);
   colz = 2.0*(omgx*vely - omgy*velx);

   rxk = relx*AX + rely*AY + relz*AZ;
   radx = relx - rxk*AX;
   rady = rely - rxk*AY;
   radz = relz - rxk*AZ;

/* centrifugal force */
   cenx = -OMEGA2*radx;
   ceny = -OMEGA2*rady;
   cenz = -OMEGA2*radz;

   source = -C_R(c,t)*( colx + cenx );
   dS[eqn] = 0;
   return source;
}

/*--------------------------------------------*/
DEFINE_SOURCE(ymom,c,t,dS,eqn)
/*--------------------------------------------*/
{
   real omgx,omgy,omgz;
   real relx,rely,relz;
   real velx,vely,velz;
   real colx,coly,colz;
   real cenx,ceny,cenz;
   real radx,rady,radz;
   real x[ND_ND],source,rxk;

   C_CENTROID(x,c,t);

   relx = x[0] - OX;
   rely = x[1] - OY;
   relz = x[2] - OZ;

   velx = C_U(c,t);
   vely = C_V(c,t);
   velz = C_W(c,t);

   omgx = OMEGA*AX;
   omgy = OMEGA*AY;
   omgz = OMEGA*AZ;

   colx = 2.0*(omgy*velz - omgz*vely);
   coly = 2.0*(omgz*velx - omgx*velz);
   colz = 2.0*(omgx*vely - omgy*velx);

   rxk = relx*AX + rely*AY + relz*AZ;
   radx = relx - rxk*AX;
   rady = rely - rxk*AY;
   radz = relz - rxk*AZ;

   cenx = -OMEGA2*radx;
   ceny = -OMEGA2*rady;
   cenz = -OMEGA2*radz;

   source = -C_R(c,t)*( coly + ceny );
   dS[eqn] = 0;
   return source;
}

/*--------------------------------------------*/
DEFINE_SOURCE(zmom,c,t,dS,eqn)
/*--------------------------------------------*/
{
   real omgx,omgy,omgz;
   real relx,rely,relz;
   real velx,vely,velz;
   real colx,coly,colz;
   real cenx,ceny,cenz;
   real radx,rady,radz;
   real x[ND_ND],source,rxk;

   C_CENTROID(x,c,t);

   relx = x[0] - OX;
   rely = x[1] - OY;
   relz = x[2] - OZ;

   velx = C_U(c,t);
   vely = C_V(c,t);
   velz = C_W(c,t);

   omgx = OMEGA*AX;
   omgy = OMEGA*AY;
   omgz = OMEGA*AZ;

   colx = 2.0*(omgy*velz - omgz*vely);
   coly = 2.0*(omgz*velx - omgx*velz);
   colz = 2.0*(omgx*vely - omgy*velx);

   rxk = relx*AX + rely*AY + relz*AZ;
   radx = relx - rxk*AX;
   rady = rely - rxk*AY;
   radz = relz - rxk*AZ;

   cenx = -OMEGA2*radx;
   ceny = -OMEGA2*rady;
   cenz = -OMEGA2*radz;

   source = -C_R(c,t)*( colz + cenz );
   dS[eqn] = 0;
   return source;
}