長岡工業高等専門学校研究紀要

論第50巻（2014）

文

仮想領域法を用いた
非圧縮性粘性流体の有限要素流体解析－ 定式化の方法の紹介と円柱周りの流れ解析 －倉橋 貴彦1・斎藤 浩一郎2・近藤 俊美3
1 機械工学科（Department of Mechanical Engineering, Nagaoka National College of Technology）
2 機械工学科（Department of Mechanical Engineering, Nagaoka National College of Technology）
3 機械工学科（Department of Mechanical Engineering, Nagaoka National College of Technology）（2014年度より：1 長岡技術科学大学

機械系）

Fluid Analysis for Incompressible Viscose Flow
Based on FEM Using Fictitious Domain Method- Introduction of Methodology of Formulation
and Fluid Analysis around Circular Cylinder Takahiko KURAHASHI1, Koichiro SAITO2 and Toshimi KONDO3
Abstract
In this study, we present some remarks on the results obtained by the FEM using the fictitious domain
method by comparing to the results obtained by the conventional FEM. As the numerical experiment,
flow analysis around circular cylinder is carried out. In the fictitious domain method, the state equation isformulated by the Lagrange multiplier method, and two meshes, i.e., foreground and background meshes,
are employed. The incompressible Navier-Stokes equation is employed as the state equation. In addition,
the fractional method is applied to calculate the state equation, and pressure and velocity fields are
independently solved. The bubble function and the linear triangular elements are applied to interpolate
velocity and pressure fields, respectively.
Key Words : fictitious domain method, finite element method, Lagrange multiplier method,
incompressible viscose flow

１．序論

一方，移動物体と流れ場のメッシュや格子を分けて計算を行う研究も行われている．有限差分法におい

流れ場中に移動物体が有る場合において流れ場の解析を行う場合，Adaptive mesh refinement法によるて用いる格子を，移動物体を副格子，流れ場を主格

1)

主格子へ反映させる重合格子法（Overset-grid法）を

移動物体近傍におけるメッシュを再生成する方法

や，回転物体の場合は，回転領域とメッシュの接合領域を分けたShear-Slip mesh update法2)が適用される．

子により表し，副格子の物理量を最小二乗法により使用した研究 3), 4)や，有限体積法に基づいた埋め込み境界法（Immersed boundary法）を使用した研究 5)

−55−

倉橋

貴彦 ・ 斎藤浩一郎 ・ 近藤

俊美

がある．また，有限要素法による計算では，仮想領域法（Fictitious domain法）が挙げられる6)- 10)．この手法ではラグランジュの未定乗数法に基づき定式化を行い，移動物体を副領域，流れ場を背景領域として表し，副領域の物理量を有限要素法において使用する補間法により背景領域へ反映させることで流れ場の解析が行われる（図-1）．仮想領域法は有限要素法に基づいた方法であるため，差分近似を主体とした方法とは異なり，背景領域の任意の点における物理量も各要素の節点の物理量から補間法により求めることができるという利点もある．以上の背景に基づき，本論文では，仮想領域法を用いた円柱周りの流れ場の解析を行う．

図-2

計算領域のイメージ図

解析モデルとしては，円柱周りの流れを対象とす
る．対象とする流体としては，非圧縮性粘性流体とし，流れ場の方程式の有限要素近似を行う要素としては，流速に対しては三角形気泡関数要素を適用し，

２．仮想領域法および有限要素法による流れ
場の方程式の定式化

圧力に対しては三角形一次要素を適用する．流れ場の方程式の解法としては，圧力場（圧力ポアソン方

計算領域のイメージ図を図-2に示す．背景領域Ω

程式）と流速場（運動方程式）を分離して解く分離に対して，式(1)，式(2)に示す非圧縮粘性流体に対型解法を使用する．有限要素方程式の数値解法としするナビエ・ストークスの運動方程式および連続式ては，圧力ポアソン方程式に対しては前処理付きのに対して仮想領域法および有限要素法を適用する．

共役勾配法，運動方程式に対しては前処理付きの安

定化双共役勾配法を適用する．本論文では，従来の有限要素法による結果との比較を通じて，仮想領域法を用いた有限要素法に基づく解析結果について考
察を行う．



1
Vi  V jVi, j  P,i 
Vi, j  V j ,i  Fi
,j
Re

(1)

Vi,i  0(2)

ここに，Viはx,y方向流速，Pは圧力，Reはレイノル
ズ数，Fiは単位体積当りの物体力を示す．初期条件
および境界条件は式(3)，(4)の様に定義する．
Vi t 0   Vˆi

(3)

in 

Vi  Vˆi on 1V x  0, V y  0 on 2 P  0 on 3(a)全体図(4)

ここに，Tiは作用力を示しており，niを外向きの単位法線ベクトルとすると式(5)のように書き表され
る．





1


Ti    P ij Vi, j  V j ,i n j
Re



in 

(5)

(b)拡大図図-1

仮想領域法におけるメッシュの一例

また，仮想領域ωにおいては，流速viαが既知である−56−

仮想領域法を用いた非圧縮性粘性流体の有限要素流体解析

ものとし，式(6)のように与える．
Vi  Viin 

 P P d* n 1
,i ,i

e

(6)



Vin,i
P*
 
d 
P,*i V jnVin, j d 
 e

t
e





流れ場の方程式を解くに際して，本研究では圧力場

と流れ場を分離して解く分離型解法を適用する．運



(10)

 P V V n d   P P n d
*

e

n n
j i, j i

n 1
,i
i

*

e

動方程式に対して発散を取り連続式を代入することにより式(7)に示す圧力ポアソン方程式が得られる．

また，式(9)の第2式の圧力項と粘性項にグリーンの定理を適用すると，式(11)のようになる．

P,iin1 

Vin,i
t



 V jnVin, j



,i

(7)

 V Vie

ここに，n，n+1は時間ステップを表す．また，運動

i

n 1



n

1

 Vi n  Vi* tV jnVi , j 2

 Vi*,i tP n 1  Vi*,j

方程式にクランク・ニコルソン法を適用すると式(8)のように書くことができる．

n  
t  n  2
Vi , j  V j ,i 2  d


Re 

1

1

e

  Vi* tFi n d   Vi * tTi n d in
e

(8)

(11)

  Vi* in 1 d

1nVin 1  Vin
 V jnVi , j 2  P,in 1

t
1
1

n  
1  n 2
n
2
V
V


j ,i   Fi
Re  i , j

, j


*

e



ここで，領域Ωeにおける積分に際して，流速に対しては三角形気泡関数要素，圧力に対しては三角形一次要素を適用し補間関数を式(12)のように定義する．
~
Vi *  N1Vi1*  N 2Vi *2  N 3Vi *3  N 4Vi 4*  N TB v *i~T
Vi  N1Vi1  N 2Vi 2  N 3Vi 3  N 4Vi 4  N B v i
 *
T *
*
*
*
 P  N1 P1  N 2 P2  N 3 P3  N p
T P  N1 P1  N 2 P2  N 3 P3  N p

ここに，n+1/2ステップの値は，n，n+1ステップの値の平均値を表す．また，副領域ωの流速の条件を考慮するために，ラグランジュの未定乗数法を導入
し，式(7)の両辺に重み関数P* ，式(8)の両辺に重み関数Vi*，式(6)の両辺に重み関数λ*乗じて，それぞれ(12)の変数の対応する領域Ωeとωeにおいて積分をすると，
重み付き残差方程式は式(9)のように得られる．
 P * P,iin1d
e

V n

  P *  i ,i  V jnVi ,nj ,i d in 


e

 t


1
 V * V n1  V n  tV nV n 2  tP n1
i
i
i
j
i
,
j
,i
e

1
 t  n 1
n  
* n 1 Vi , j 2  V j ,i 2  d  e Vi i dRe, j n*  Vi t Fi d in  e *i Vi n1  Vi d  0 in  e

点の番号，4は要素内重心点の番号を示す．面積座標を(η1 ，η2 ，η3)とすると，N1=η1 ，N2=η2 ，N3=η3 ，
N4=27η1η2η3と表される．また，領域ωeにおける積分







ここに，Niは形状関数を示し，添え字の1,2,3は各頂



に際して，補間関数を式(13)のように定義する．
(9)

*i  *i x  x  a i  i x  x  a (13)

ここにδはデルタ関数を示している．式(13)を使用



して，式(9)，(11)における領域ωeの積分を行うと式(14)，(15)のように書くことができる．

ここにλiはラグランジュの未定乗数を示す．式(9)の二階微分の項にグリーンの定理を適用すると，式(10)のようになる．

−57−



e

Vi* i d NdV a a
i 1

*
i

i

(14)

倉橋



e

*i Vi  Vi d Nd

*i aVi a 

i 1

Nd



貴彦 ・ 斎藤浩一郎 ・ 近藤

*i aVi a俊美 V V  V V



n n
n n
n 1
nx
x x, xy x , y  P, xe
 V ynV xn, y  V ynV yn, y  P, ny1 n y Nd

b

(15)



i 1

p  P1

P3 

P2


V x3
v y  V y1 V y 2 V y 3

ここに，Ndは対象とする要素Ωe内に存在する副領域v x  V x1 V x 2の節点の数を示し，aは，副領域の節点座標を示す．
ここで，Vi*(a)，Vi(a)を要素Ωe 内の補間関数により

 

T

(18)




~ T
Vx4
~ T
V y4

表すと，式(16)のようになる．
また，図-3のように要素領域Ωe内に副領域ωeの節点
~*
*
*
*
*

Vi a   N1a Vi1  N 2 a Vi 2  N3 a Vi3  N 4 a Vi 4~
Vi a   N1a Vi1  N 2 a Vi 2  N3 a Vi3  N 4 a Vi 4(16)

今，図-3のように，対象とする要素Ωe内に副領域の節点が1点ある状態を想定し，有限要素方程式を導く．ここに，副領域の節点座標はa=(xa,ya)とする．
まず，圧力ポアソン方程式は副領域の値には直接関係がないため，従来からの有限要素法に基づく定式化と同様に式(17)に示すようになる．

C  C p
xx

yy

n 1

 1
 
S x v nx  S y v ny  D xx  D xy v nx
 t D yx  D yy v ny  b n





(17)



が存在する場合，運動方程式に対する有限要素方程
式は式(19)に示すようになる．
n 1  t
H xx

M 0 0 0   v x  v 
 0 M 0 0  y 
 0 0 0 0  r x  0 0 0 0  r y M 0 0 0   v x  
0 M 0 0  v y  0 0 0 0  r x  0 0 0 0  r y 
tA x p tA y p 


 0 
 0 

n 1

2
 t H  2 yx
 Bx
 0


n 

 t H xx 2tH  2 yx
0

 0

t H
xy
2
t H
yy
2

B Tx

0

0

By

0

t H
xy
2
t H
yy
2

0

0

n 

 tt x tt y  fx 
 fy 00 v x  
B Ty   v y  
0  rx 
0   r y 

0 0  v x  
0 0  v y  
0 0  r x 0 0  r y 

n 

(19)
ここで，式(17)の係数行列は式(18)のように書くことができる．

ここに，式(19)における各係数行列は式(20)のように書くことができる．

 N N dC 
 N N dS  NN dS  NN dD 
 N N v N dD 
 N N v N dD 
 N N v N dD 
 N N v N d
C xx 
yy

x

y

xx

xy

yx

yy

e

e

e

e

e

e

e

e

,x

T
,x

,y

T
,y



M  N B N TB d
e

H xx

 N N v N d N N v N d
N N d  2
 N N d  H  N N dH 
 N N dH 
 N N v N d N N v N d 
N N d  2
N N d 





T
B, x



T
B, y

,x

,x

,y

,y

T
B

T
B

T
B

T
B

x

y

x

y

e

e

T
B, x

1
Re

T
B, y
T
B, x

T
B, y

T
B

x

T
B, x

B

T
B

y

T
B, y

B, x

T
B, x

e

xy

1
Re

yx

1
Re

yy

e

1
Re

−58−

B

e

B

e

e

e

B, y

T
B, x

B, x

T
B, y

T
B

x

B

T
B

B, x

y

e

B, y

T
B, y

B, y

T
B, y

T
B, x

T
B, y

T
B, x

e

n 1

仮想領域法を用いた非圧縮性粘性流体の有限要素流体解析

x3, y3 

 N N dA 
 N N d

Ax 

B, x

e

y

B, y

e

T

T

B x  N 1 a  N 2 a  N 3 a  N 4 a x1, y1 

B y  N 1 a  N 2 a  N 3 a  N 4 a 

p  P1




P2

P3 T




~ T
v x  V x1 V x 2 V x3 V x 4~ Tv y  V y1 V y 2 V y 3 V y 4

r x   x a , r y   y a ea  xa , ya x2 , y2 

 T Nd, t   T Nd

図-3

ここに副領域に働く流体力tx ，ty は求まった流速，

３．数値実験

tx 

e

x

y

e

e

(20)

y

f x  V x a , f y  V y a 

背景領域と副領域の重合部分の一例

圧力の値を用いて式(21)により計算することができ
る．

３．１

レイノルズ数 250，時間増分量 0.001 の場合に対する数値実験

図-4に示す解析モデルに対して，円柱周り流れのt x n  1 M v x n 1  v x n tt H xx  v x n 1  v x n2t
H xy v y n 1  v y n  t A x pn 1 

2

1
n
n 1
n
M  v y  v y
ty 
t
t
H yx v x n 1  v x n2tH yy v y n 1  v y n  t A y pn 1 

2


 



                  

解析を行う．本論文では，従来の有限要素解析結果と，仮想領域法を用いた有限要素解析結果の比較を行う．各解析においては表-1,2に示す計算条件を用
(21)

いる．
t x=0, v y=0(vx,vy)=(1.0,0.0)

1.0

Fictitious dimain FEM(vx,vy)=(0.0,0.0) in domain
Conventional FEM(vx,vy)=(0.0,0.0) on boundary

また，背景領域の各要素における副領域の各節点のt x=0, v y=0

要素内に副領域の節点（節点座標((xa,ya)）が存在する場合を考えると，式(22)に示す行列式D1, D2, D3の副領域の節点があることを判定できる．

x1  x a D1 
x2  xa
x2  xa D2 
x3  x a D  x3  x a
 3 x1  x ay2  yay3  y a

図-4

時間増分量Δt
節点数
要素数

(22)

y1  y a

背景領域 ω の節点数
背景領域 ω の要素数
背景領域 Ω の節点数
背景領域 Ω の要素数

−59−

250
0.001
1,636
3,116

表-2 仮想領域法を用いた有限要素解析の計算条件
レイノルズ数
時間増分量Δty3  y a

計算モデル図表-1 従来の有限要素解析の計算条件
レイノルズ数y1  y ay2  ya

8.0

Center of circle(X,Y) = (0,0)

探索を行う場合は，図-3のように背景領域のある要
素の三点の座標を(x1,y1)，(x2,y2)，(x3,y3)とし，その値が全て正の値となり，背景領域内の要素において

p=0

12.0

250
0.001
331
600

1907
3716

倉橋

貴彦 ・ 斎藤浩一郎 ・ 近藤俊美

解析結果を図-5～8に示す．図-5,6はT=50における仮想領域法を用いた有限要素法により得られた圧力
分布・流速ベクトル図と，従来の有限要素法により得られた圧力分布・流速ベクトル図を示す．(a)は解析領域全体の結果， (b)は円柱周りの結果の拡大図を示す．両図の比較より，同等な流れ場が得られていることを確認できる．また，図-7,8はそれぞれ
計 測 点 (x,y)=(4.98,-0.98) に お け る 流 速 U（U=(u2+v2)1/2）の値，圧力Pの値の時系列変化を示す．

実線は仮想領域法を用いた有限要素法による結果，
点線は従来の有限要素法による結果を示す．流速U
および圧力Pの値とも，位相が若干ずれてはいるものの同等な結果が得られていることがわかる．
(a)解析モデル全体図(a)解析モデル全体図(b)円柱周りの拡大図
図-6T=50における圧力分布図と流速ベクトル図（従来の有限要素法(Re=250，Δt=0.001)）

図-7(b)円柱周りの拡大図
図-5T=50における圧力分布図と流速ベクトル図

測定点（x,y）=(4.89,-0.98)における流速U（U=(u2+v2)1/2）の値の時系列変化（仮想領域法を用いた有限要素法(Re=250，Δt=0.001)）−60−(Re=250，Δt=0.001)

仮想領域法を用いた非圧縮性粘性流体の有限要素流体解析

図-8

３．２

測定点（x,y）=(4.89,-0.98)における図-10

測定点（x,y）=(4.89,-0.98)における圧力Pの値の時系列変化

圧力Pの値の時系列変化(Re=250，Δt=0.001)(Re=200，Δt=0.001)

レイノルズ数 200，100，時間増分量 0.001の場合に対する数値実験

次にレイノルズ数を200と100に変えた場合の結果を示す．計算条件は前節と同様の条件とする．図-9～12はレイノルズ数を200と100の場合に対する計測点(x,y)=(4.98,-0.98)における流速 U （ U=(u2+v2)1/2 ） の値，圧力Pの値の時系列変化を示す．前節において示したレイノルズ数が250の場合に比べ，レイノルズ数が小さくなるにつれて，仮想領域法を用いた場
合の有限要素法による結果と従来の有限要素法の結果の位相差は小さくなっていることがわかる．このため，レイノルズ数が高くなるに従い，副領域の情図-11

測定点（x,y）=(4.89,-0.98)における報を背景領域に伝える際の補間誤差を軽減できるよ

流速U（U=(u2+v2)1/2）の値の時系列変化

うに，よりメッシュを細かくする等の必要があると(Re=100，Δt=0.001)

考えられる．

図-12
図-9

測定点（x,y）=(4.89,-0.98)における測定点（x,y）=(4.89,-0.98)における圧力Pの値の時系列変化

流速U（U=(u2+v2)1/2）の値の時系列変化(Re=100，Δt=0.001)(Re=200，Δt=0.001)

−61−

倉橋

３．３

貴彦 ・ 斎藤浩一郎 ・ 近藤

俊美

レイノルズ数 250，時間増分量 0.01，
0.0001 の場合に対する数値実験

次に3.1節の計算条件のもと，時間増分量を0.01，
0.0001に変えた場合の結果を示す．図-13～16は計測点(x,y)=(4.98,-0.98)における流速 U （ U=(u2+v2)1/2 ）の値，圧力Pの値の時系列変化を示す．3.1節の時間
増分量の条件（Δt=0.001）に比べて時間増分量の設定を大きく（Δt=0.01）した場合は，図-7,8の場合に比べて位相差が生じていることを確認できる．また，図-7,8の場合に比べて，時間増分量を小さく（Δt=0.0001）設定した場合は，図-7,8に示す結果と同様に位相差を完全に無くすことはできない結果

図-15

測定点（x,y）=(4.89,-0.98)における流速U（U=(u2+v2)1/2）の値の時系列変化となった．

(Re=250，Δt=0.0001)

図-13

測定点（x,y）=(4.89,-0.98)における流速U（U=(u2+v2)1/2）の値の時系列変化

図-16(Re=250，Δt=0.01)

測定点（x,y）=(4.89,-0.98)における圧力Pの値の時系列変化(Re=250，Δt=0.0001)

４．結論
本研究では，非圧縮性粘性流体に対する円柱周り
流れの解析に仮想領域法を用いた有限要素法を適用し，従来の有限要素法による結果との比較を行った．
得られた知見を以下に整理する．
・流速ベクトルおよび圧力分布の図より，仮想領域法を使用した有限要素法と従来の有限要素法の結果
図-14

測定点（x,y）=(4.89,-0.98)におけるには大きな差がないことがわかった．

圧力Pの値の時系列変化(Re=250，Δt=0.01)

・ある測定点における流速U（U=(u2+v2)1/2）と圧力Pの値について，仮想領域法を使用した有限要素法と従来の有限要素法の結果の比較をしたところ，レイノルズ数が大きくなるに従い，従来の有限要素法の結果と仮想領域法を用いた有限要素法の結果に位−62−

仮想領域法を用いた非圧縮性粘性流体の有限要素流体解析

Method,14th

相差が生じることを確認できた．

Symposiumon

Coputational

Fluid

Dynamics,C01-2，(2000), pp.1-5.

・レイノルズ数が250，時間増分量0.001の条件にお4)method for 3D unsteady incompressible flows, J. Comput.

いて，時間増分量Δtを10倍（Δt=0.01）に変えたと
ころ，位相差が大きくなる結果が得られた．逆に，
時間増分量Δtを1/10倍（Δt=0.0001）した場合は，Tang,H.S., Jones,S.C., Sotiropoulos,F., An overset-gridPhy., Vol.191, (2003), pp.567-600.

5)Kino,C., Watanabe,T., Nishida,A. and Takemiya,H., Flowanalysis around an in-line forced oscillating circular

時間増分量が0.001の場合と然程変わらない結果とcylinder using IB-method, Transactions of the Japanなった．

Society of Mechanical Engineers Series B,Vol.78, No.796,(2012), pp.60-73.

謝辞：本論文を執筆するにあたり長岡技術科学大学

6)

Glowinski,R., Pan,T.W., Hesla,T.I., Joseph,D.D., andPeriaux,J., A Fictitious Domain Method with Distributedより「H25年度 高専－長岡技科大共同研究助成」のLagrange Multipliers for the Numerical Simulation of

援助を受けた．ここに謝意を表す．また，本論文のParticulate Flow, Contemp. Math., Vol.218, (1998),

解析結果は，九州大学情報基盤研究開発センターのpp.121-137.

高性能アプリケーションサーバSR16000を使用し計

7)

Glowinski,R., Pan,T.W., Hesla,T.I., Joseph,D.D. and算を行ったものである．計算機の使用についても，

Periaux,J., A Fictitious Domain Approach to the Direct

センター関係者の方へ謝意を表す．本研究は，高Numerical Simulation of Incompressible Viscous Flow

専・両技大間教員交流制度による長岡技術科学大past Moving Rigid Bodies: Application to Particulate Flow,

学・長岡工業高等専門学校間における人事交流中の内容であり，長岡工業高等専門学校の機械工学科のCommun. Comput. Phys., Vol.169, (2001), pp.363-426.
8)with Distributed Lagrange Multiplier for Particulate Flow,

教職員からは，研究成果発表会等を通じてたくさん

Int. J. Comp. Meth. in Engi. Sci. and Mech., Vol.8, (2007),の建設的な意見を頂いた．合わせてここに謝意を表す．

pp.115-122.
9)

Vanella,M.,

Rabenold,P.,

embedded-boundary

Balaras,E.,A

methodElastodynamics, SIAM J. SCI. COMPUT., Vol.29, No.3,

direct-forcing

withadaptivemeshrefinement for fluid-structure interaction problems, J.

(2007), pp.1244-1267.
10) Glowinski,R., and He,Q., A Least-Squares/Fictitious

Comput. Phy., Vol.229, pp.6427-6449, 2010.
2)

Behr,M. and Tezduyar,T.,The Shear-Slip Mesh Update
Method,Comput. Meth. Appl. Mech. and Eng., Vol.174,(1999), pp.261-274.

3)

Circular

Domain Method for Linear Elliptic Problems with Robin
Boundary Conditions, Commun. Comput. Phys., Vol.9,
No.3, (2011), pp.587-606.
(2014. 3. 31 受付)

Ogasawara,K. and Kuroda,S.,Flow Analysis Around
Oscillating

Becache,E., Rodriguez,J., and Tsogka,C., A Fictitious
Domain Method with Mixed Finite Elements for

参考文献
1)

Nagai,M., and Kawahara,M., A Fictitious Domain Method

CylinderUsing

OversetGrid−63−