日本機械学会［No.167-1］北陸信越支部 第 53 期総会・講演会 講演論文集 [2016.3.5 長野県長野市]
［No.167-1］日本機械学会北陸信越支部 第 53 期総会・講演会 講演論文集 [2016.3.5 長野県長野市]
306
仮想領域有限要素法を用いた物体周りの流れ解析による流体力の評価
寺門 幸宏*1，倉橋 貴彦*2
Evaluation of fluid force based on flow analysis around body using fictitious domain FEM
Yukihiro TERAKADO
*1and Takahiko KURAHASHI
*2*1 Graduate school of Nagaoka University of Technology, 1603-1 Kamitomiokamachi, Nagaoka-shi, Niigata, Japan
When the fluid analysis is carried out in the case that there is a moving object in the flow field, the mesh regeneration method is applied due to the adaptive mesh refinement method. In addition, the shear-slip mesh update method is employed in case of rotating body problems. On the other hand, a method using the flow domain overlapped with the moving object is also proposed. In the finite difference method, a sub-grid is represented by moving object, and a main- grid is used for flow field calculation. The overset-grid is applied to the fluid analysis to reflect the sub-grid to main-grid by the least-squares method. Furthermore, the immersed boundary method based on the finite volume method is also employed. In the finite element analysis, the fictitious domain method is often adopted for the moving body problems. In this method, the computational domain is divided into foreground domain and background domain. In addition, the formulation is carried out based on the Lagrange multiplier method, and is applied to consider the velocity condition in the foreground domain. The finite element fluid analysis is carried out to reflect the physical quantity of the foreground domain to background domain. Physical quantities at any points of the foreground domain can be obtained from the physical quantities at each node by interpolation method. In this study, the flow field analysis using the fictitious domain finite element method is carried out.
Key Words : Fluid analysis, Finite element method, Fictitious domain method, Lagrange multiplier method, Interpolation method
1. 緒 言
有限要素法を適用した解析手法は，航空機やロケット等の流れ解析や強度計算などにも用いられ，汎用性の高 さから様々な分野で利用されている．しかし形状最適化問題(1)(2)など流れ場内の物体を変形させながら解析を行 う場合には，その物体形状に合わせてメッシュを再構築する必要があるため，それに伴う計算コストの増加が問 題となっている．そこで，計算領域を流れ場の背景領域と移動する物体の前景領域に分離した解析手法の研究が 注目されている．仮想領域有限要素法(3)(4)は領域を分離した解析手法のひとつであり，前景領域の持つ物理量をラ グランジュの未定乗数法により背景領域内に作用させて解析を行う手法である．ラグランジュ未定乗数法は制約 条件付きの関数の極値を求める数学的な手法であり，流れ場の支配方程式に対する付帯条件として計算すること ができる．支配方程式はナビエ・ストークスの運動方程式と連続式を適用する．分離型解法より導いた圧力ポア ソン方程式と流速場における運動方程式を離散化し，本手法における補間関数を考慮することで前景領域の流速 の条件を含んだ式が得られる．本研究では非圧縮粘性流体について流れ場内の物体に及ぼす流体力を解析し，従 来法との比較によって本手法の実用性を検証する．また，前景領域の形状を変更した際の物体周りにおける流体 力の変化を示す．
*1 長岡技術科学大学大学院 機械創造工学専攻（〒940-2188 新潟県長岡市上富岡町 1603-1）
*2 正員，長岡技術科学大学大学院 機械創造工学専攻（〒940-2188 新潟県長岡市上富岡町 1603-1）
E-mail: [email protected]
2. 弱形式の離散化 2・1 非圧縮粘性流体に対する流れ場の支配方程式
背景領域Ωの流体における流れ場の支配方程式として式（1）にナビエ・ストークスの運動方程式，式（2）に 連続式を示す．これらの式に対して仮想領域法を用いた有限要素法を適用し，流れ場の解析を行う．ここで，uiは 流速，Pは圧力，Reはレイノルズ数，fiは単位面積あたりの物体力を示す．
ij ji
j i ij i j
i uu P u u f
u
, , , , ,
1
（1）
,i 0
ui （2）
初期条件と境界条件を以下に示す．前景領域ωにおいては流速が既知であり、式（5）で定義する．
0 ˆi in i t u
u （3）
3 2 1
on 0
on 0 , 0 ˆ on P
u u
u u
y x
i i
（4）
in
i
i u
u （5）
2・2 有限要素方程式の導出
本研究では分離型解法を適用し、圧力場と流れ場を分離して解く．式（1）についてオイラー法を使用し時間方 向の離散化を行うと式（6）となる．
j inn i j n
j i n i n
j i n j n i n
i u u P u u f
t u
u
, , , 1 , ,
1 1
（6）
式（6）を微分し，発散を取って連続式を代入し定数項と高次の微分項を消去することで，以下に示す圧力ポアソ ン方程式が得られる．
i n
j i n j n
i i n
ii u u u
P t
, , ,
1 ,
1 （7）
次に，式（1）の運動方程式に対してクランク・ニコルソン法を適用し，流れ場についての離散化を行う．両辺 に重み関数u*iを乗じ，背景領域に対して積分をすることで式（8）に示す重み付き残差方程式が得られる．n，n+1 は時間ステップであり，n+1/2はn，n+1の平均値を表す．
e
n i e i
j n
i j n
j i n
i n
j i n j n i n i
i tP t u u d u tf d
u tu u u
u *
, 2 1 , 2 1 , 1
, 2 1 , 1
*
（8）
ここに，圧力に対しては三角形一次要素，流速に対しては三角形気泡関数要素を補間関数として適用すること で，空間方向に離散化した式を得ることができ，有限要素方程式（式（9））を得る．ここでCは粘性行列，S は移流行列，Mは質量行列，Hは移流行列と粘性行列の足し合わせたものを表す．
1 in
1 1
n i n i ij n
i ij
n n i ij i n ii
p A u H M u
H M
b u D S p
C （9）
2・3 前景領域を考慮した離散化
ここで，背景領域内に物体が存在する場合を考える．Glowinski らの研究(4)に基づき，式（5）の制約条件を考 慮し，式（9）の第2式を解くことを考えると，式（10）第1式に示す有限要素方程式が得られる．また式（10）
の第2式に示すaは前景領域の節点座標を表す．
a u u B
p A u H M r
B u H M
i n i i
n i n i ij n
i T i n i ij
ˆ
in
1
1 1
1
（10）
ここに式（10）の第2式は有限要素法における補間関数を示しており，右辺のハット付き変数は既知量であるこ とを示す．また，ラグランジュの未定乗数は従来の有限要素法におけるトラクションベクトルを表す(5)．このこ とより，ラグランジュの未定乗数を足し合わせることで物体に働く流体力を簡便に導くことができる．
3. 仮想領域法を用いた数値実験
3・1 円柱周りの解析における仮想領域有限要素法および従来の有限要素法による結果の比較
円柱周りの流れ場に対して，従来の有限要素法を使用した解析と仮想領域法を用いた有限要素解析の比較を行 う．解析モデルと計算条件をそれぞれ図1，表1に示す．本実験では密度等を考慮した次元（単位）を考慮した 解析を行う．従来法での流体力TAの計算は式（11），（12）に示すように，1ステップ間での流速ベクトルの差を 計算し，物体周りの要素（総数:Na）における値を足し合わせることにより求めることができる．一方，仮想領域 法の解析による流体力 TBは式（10）よりrin1と等価になるため，前景領域と重なる要素（総数:Nb）におけるラ グランジュ未定乗数の足し合わせによって計算することができる（式（13））．
1
1
ij ni ij ni i n
n
i M H u M H u Ap
t （11）
a
i N
k k n i n
TA 1
t （12）
b i
N
m m n i n
TB 1
r 1 （13）
以下の図2，3に従来法と仮想領域有限要素法におけるT=10 [s]での圧力分布と流速ベクトルを示す．仮想領域有 限要素法の場合には物体の内部にも要素が存在するが，制約条件より前景領域の流速が考慮されるため，流速ベ クトルは従来法と同様に円柱周りに分布する．また，従来法と本手法について抗力の時系列変化を図4に示す．
Density [kg/m3] 1.00
Kinematic viscosity [m2/s] 0.01 Time increment [sec] 0.001 Sectional area of domain ω [m2] 0.7854
Table 1 Computational conditions
Fig.1 Computational model
3・2 前景領域の形状を変えた場合に対する流体力の変化
仮想領域有限要素法では二つの領域を合わせて計算するため，条件に合わせて内部の物体を自由に変更するこ とができる．この特徴を利用して，流体中の物体の形状だけを変える処理を付加した．それに伴う各節点の移動 量を求めるにあたり，式（14）に示すラプラス方程式に対して有限要素法を適用する．境界条件には式（15）の ディリクレ条件を与える．また本実験の制約条件として物体の断面積を一定とした．
0 0
2 2 2 2
2 2 2 2
y y x
y y
x x
x
（14）
in
ˆ
ˆ
i new i
i new i
y y
x
x （15）
いま，流れ場内の物体を変更するため節点の移動について考える．支配方程式と境界条件より前景領域の各節点 の移動量Δx，Δyを計算し，Δx，Δyの値を元の座標値x，yにそれぞれ足し合わせることで節点を移動させる．そ の後，再構成された前景領域の形状を用いて流体解析を行う．物体は楕円率（短軸半径／長軸半径）が0.64（Ellipse1）
と0.44（Ellipse2）の2つの楕円形状を使用し，抗力の違いについて検証を行った（図5）．楕円形状と円柱の抗力
を比較した結果を，図6に示す．結果から分かるように，円形から楕円率が小さい楕円になるにつれて抗力は減 少した．解析終端時刻において楕円形状の抗力を円柱の場合と比較したところ，Ellipse 1は18%，Ellipse 2では 30%低下した．
Fig.2 Pressure and velocity vector distributions in case of conventional FEM.
Fig.3 Pressure and velocity vector distributions in case of fictitious domain FEM.
Fig.4 Drag force around the cylinder comparing with a conventional FEM and the Fictitious domain FEM.
4. 結 語
本研究では，仮想領域有限要素法を適用した流体解析によって物体に作用する流体力を評価した．本手法によ り背景領域を変更することなく様々な物体について解析を行えることが確認できた．計算精度においても流体力 の比較から従来の有限要素法と同等の結果を得られたため，本手法の妥当性を実証することができた．また本研 究ではラプラス方程式を利用した形状の変更と，その物体を使った流れ解析により抗力の違いについて検証を行 った．今後の展望として，形状最適化の手法を本研究に取り入れ，翼形状の最適設計を目指す．ラプラス方程式 による有限要素法を利用して形状を更新し，その物体に対して仮想領域法による流体解析を施す．これらの手法 を導入することで，前景領域の流体解析を繰り返すような最適形状問題において効率的な解析が期待できる．
謝 辞
本研究の実施にあたり，金沢工業大学空力設計研究室の佐々木大輔講師，並びに長岡技術科学大学機械創造工 学専攻航空流体工学研究室の山崎渉准教授に流体中における形状最適化に関する資料等を頂いた．ここに記して 謝意を表す．
文 献
(1) Yamazaki, W., “Advanced Kriging response surface approach for efficient design optimization”, Transactions of the JSME, Vol.
80, No. 818 (2014), p. TRANS0285.
(2) Sasaki, D., Takeguchi, Y., Obayashi, S., Nakahashi, K., and Hirose, N., “Multi-Objective Aerodynamic Optimization of the Next Generation SST Wing for Supersonic and Transonic Cruise”, Journal of the Japan Society for Aeronautical and Space Sciences, Vol. 47, No. 551 (1991), pp. 464-469.
(3) Kurahashi, T., Saito, K., and Konda, T., “Fluid Analysis for Incompressible Viscose Flow Based on FEM Using Fictitious Domain Method -Introduction of Methodology of Formulation and Fluid Analysis around Circular Cylinder-”, Research reports of the Nagaoka Technical College, Vol. 50 (2014), pp. 55-63.
(4) Glowinski, R., Pan, T.W., Periaux, J., “A fictitious domain method for Dirichlet problem and applications”, Computer Methods in Applied Mechanics and Engineering, Vol. 111 (1994), pp. 283-303.
(5) Glowinski, R., Pan, T.W., Kearsley, A.J., and Periaux, J., “Numerical Simulation and Optimal Shape for Viscous Flow by a Fictitious Domain Method”, International Journal for Numerical Methods in Fluids, Vol. 20 (1995), pp. 695-711.
Fig.5 Cylinder model and 2 types of Ellipses mesh model on the background mesh.
Fig.6 Drag force around the cylinder and Ellipses comparing with a conventional FEM and the Fictitious domain FEM.
