第 27 回数値流体力学シンポジウム
C02-3

Immersed boundary 法における簡便な流体力算出方法について
A simple hydrodynamic force evaluation method for immersed boundary methods⃝

野々村 拓，宇宙機構宇宙研, 神奈川県相模原市, E-mail：nonomura@ﬂab.isas.jaxa.jp
大西 順也，東大生研,兵庫県神戸市, E-mail：jonishi@iis.u-tokyo.ac.jp

藤井 孝藏，宇宙機構宇宙研, 神奈川県相模原市, E-mail：nonomura@ﬂab.isas.jaxa.jp

Taku Nonomura, ISAS/JAXA, Sagamihara, Kanagawa, Japan.
Jyunya Onishi, IIS, University of Tokyo, Kobe, Hyogo, Japan.
Kozo Fujii, ISAS/JAXA, Sagamihara, Kanagawa, Japan.
The simple hydrodynamic force evaluation method for immersed boundary method is proposed and examined.
The analytical solution of the Stokes ﬂow is adopted and the error in the ﬂow ﬁeld itself is perfectly eliminatedto purely discuss the force evaluation method. The proposed method which is based on simple integration on
the grid face works well compared with the polygon-base method. This is because the surface area vector of
the grid surface can be converged to the analytical surface area vector of the object which is modeled in theimmersed boundary method. Also, the appropriate evaluation of stresses on the grid face is required for theaccurate estimation. This can be achieved by the simple second order diﬀerence scheme using the cell inside theobjects.

諸言Immersed boundary 法は複雑形状周りや混相状態の非圧縮／圧縮流体の解析方法として，注目を集めている．
(1, 2, 3, 4, 5)
本論文では，Immersed boundary 法における流体力評価に関して論ずる．
一般的に流体の数値解析において物体にかかる流体力の計算は以下のように行われる．
∫Fj =σij ni dS,
(1)

1.

ΓS

ここで σij は応力テンソル, nj は面に垂直な単位ベクトルである．応力テンソルは，以下のように定義される：()∂uj∂uiσij = −pδij + µ+,
(2)
∂xi
∂xj

これらは，曲線座標を用いた場合は直接的に行われているが，直交格子を利用した Immersed Boundary 法など
では，物体は格子で表現されておらず (Fig. 1) ，その一般的な評価方法は未だ確立されていない．そこで，本発
表では，その評価方法に関して議論を行った．特に，格子面だけの積分だけで流体力を算出するシンプルなモデルを提案し，その精度を論ずる．
これまでの流体力の評価方法の精度に関しては，あまり議論されておらず，数値解析の結果と同時に評価されることが多い．その際の誤差には，数値解析の結果に於ける誤差と流体力の評価方法における誤差が含まれるため，流体力の評価方法のみから生ずる誤差を取り出して注意深く論ずる必要がある．本論文では，厳密解が存在する Stokes 流れを利用することで，数値解析の結果に含まれる誤差を完全に取り除き，流体力の評価方法の誤差のみを取り出して議論することを狙った．
計算方法
本節では，2.1 において計算格子の説明およびセルの定義を行う．その後，2.2 で一般的に行われていると考えられるポリゴンベースの流体力評価方法，2.3 で本論文で提案するメッシュベースの流体力の評価方法を提案する．

2.

ここで p は圧力, ui は速度の各方向成分， µ は粘性係数を示す．このような評価方法は Near-ﬁeld 法(6) と呼ばれ，圧力抵抗および粘性抵抗を分解を行った上で非定常力を正確に算出するのに使われる．一方, 検査体積を遠方場に設ける Far-ﬁeld 法や空間を積分する Mid-ﬁeld 法もある．(6) Far-ﬁeld 法は，物体を取り囲む仮想的な領域を考え，その領域に出入りする運動量の変化から，物体に作用する流体力を評価する手法である．Immersedboundary 法に対する Far-ﬁeld 法の適用例としては，Laiと Peskin(7) により提案された手法（control volume 法）があるが，この手法は定常状態にのみ適用が可能である．
その後，Balaras(8) によって非定常問題，Shen ら(9) によって移動物体問題への拡張が行われた．ただし，Far-ﬁeld法や Mid-ﬁeld 法で求められる流体力（特に抵抗力）は，造波抵抗，形状抵抗，誘導抵抗などに分解され，単に圧力成分と粘性成分に分けることができない．したがって，流体力の総和を解析することが目的であれば，どちらの手法を用いてもよいが，そうでないのであれば，解析の目的に応じた使い分けが必要である．これらの方法は非定常力の算出が困難で，計算コストの面でも負荷が大きいためここでは考慮しない．
Eq. (1) を用いて，数値解析で力を求めるには以下の 2点が必要となる：• 表面の離散化• 表面での物理量の評価
1

2.1 セルの分類
計算格子を Fig. 2 に示す．図にあるように流体領域および固体領域を含む全領域で等間隔直交格子を用いており，格子によって分割されたセルは，物体境界との位置関係で以下の 3 グループに分けられる：• 流体セル• 固体セル• 境界セル .
本分類は，物体形状が複雑な場合には特に，必ずしも自明でないが，この点に関しては本発表の目的の範疇を超えるため，ここでは議論を行わない．また固体セルの内，後述する流体計算で物体内部の流体量を仮想的に定義する点をゴーストセルとする．
流体力を Eq. 1 で計算するために必要な圧力，速度などの流体情報は，流体セルで定義されており，物体内部(ゴーストセル) で流体量が必要な場合は内挿や外挿を行う必要がある．本発表では，セル中心での定式化を論ずるが，セル頂点での定式化への適用は容易である．
2.2 ポリゴンベース法
本方法 (以下 PM： Polygon-based method) は，流体力の評価方法として一般的と考えられるものである．本方c 2013 by JSFM
Copyright ⃝第 27 回数値流体力学シンポジウム
C02-3

ﬂuid%cell
boundary%cell
solid%cell

ﬂuid%regionsolid%region

Fig. 2: The classiﬁcation of the computational cells.

法では，物体形状は 3 次元のケースではポリゴンの集合として定義，2 次元のケースでは線分の集合として定義される．レベルセット法等で界面が定義されている場合も，セルと界面の交線が精度に応じて計算でき，物体形状を表すポリゴンを作成することができる．これ以降ポリゴン，線分といったこれらの表面要素を示すのに “facet” を用いる．
本方法は，直接的に表面積分 Eq. (1) を，facet を用いて離散化する：∑Fj =(σij ni ∆S)facet ,(3)facet(a) Body-conforming grid

ここで ∆S は facet の表面積であり，積分は全 facet に対して行う．この積分を 2 次精度で行うためには中点則を用いればよく，括弧内の量を facet の面中心で求めれば良い．括弧内は，幾何学的な面積，単位垂線ベクトル，と流体量である応力に分けられる．ここで，前者に関しては，facet の頂点情報から求めることができる．一方流体の応力は，の方法に準じて後述するイメージポイントを利用して求める．
本方法の模式図を Fig. 3(a) に示す. facet の中点から垂線を伸ばしイメージポイント I1 ， I2 をこの垂線上の距
離 d1 および d2 の位置に作成する．これらのイメージポ
イントは，facet 中点での応力を求めるのに利用される．
圧力に関しては，線形補間を用いる．

p(d) = C1 d + C0 .

(b) Non-conforming grid(4)

ここで，係数 C1 と C0 は以下の条件を満たすように決定する：
p(d1 ) = p1 , p(d2 ) = p2 ,(5)

Fig. 1: Body-conforming and non-conforming grids.

ここで p1 と p2 はそれぞれ，参照点 I1 と I2 における圧力値である．結果として，以下の表式が得られる：

pI = p(0) = p1 −d1(p2 − p1 ).
d2 − d1(6)次にせん断応力に関しては速度の勾配が必要であり，これを求めるために速度の補間が必要となる．ここでは，放物線補間を用いる．

u(d) = C2 d2 + C1 d + C0 ,2(7)c 2013 by JSFM
Copyright ⃝第 27 回数値流体力学シンポジウム
C02-3

ここで，係数 C2 ,C1 ,C0 は以下の条件を満たすように決定する：

u(d1 ) = u1 ,u(d2 ) = u2 ,u(0) = uB ,(8)ここで u1 と u2 はそれぞれ，参照点 I1 と I2 のおける速度であり, uB は facet 中点での速度である．ここで，facet中点での速度は滑り無し，透過無し条件のもとで物体の速度の速度に等しい．この条件のもと，点 I での速度勾配が以下のように求まる：

u1 d22 − u2 d21
∂u=,
∂n I
d1 d2 (d2 − d1 )∂u= 0,∂s I

probe

I2∂u= 0, (9)∂t I

ここで座標 n, s および t はそれぞれ facet に垂直，平行方向を示す．他の速度成分に関しても同様に計算できる．.
上述の内挿操作を行うためには，流体量をイメージポイントに補間する必要がある．本発表では，この補間には 3 重 (2 次元の場合 2 重) 線形補間を用いた．ここで，イメージポイントまでの距離を再帰参照を防ぐために d1と d2 はそれぞれ 1.75∆x と 3.5∆x に設定した．

2.3 格子ベース法
上記のポリゴンベース法に替わる方法として格子を用いた簡便な方法 (以下 MM： Mesh-based method) を提案する．ここでは，本技術のキーをわかりやすく示すため，適切な評価をしたものや，不適切な評価をしたものの 2 手法を取り上げる．これら 2 手法は積分の方法は同じであるが，用いる応力の評価方法が異なる．まずは，2手法（後述する MM1,MM2）に対して共通の積分方法を示す．模式図を Fig. 3(b) に示す．本方法では Eq. (1) の面積分を各方向成分に分解し, それらを対応する格子面を用いて近似する．
∫∫
∫
Fx =
σxx nx dS + σyx ny dS + σzx nz dS (10)∑=(σxx nx ∆S)face,xface,x∑+(σyx ny ∆S)face,yface,y∑(σzx nz ∆S)face,z .
(11)+face,z= pi,j ,= pi,j ,d1

Id2

facet&center
points&for&interpola.on

(a) Polygon-based method

face,xIxここで，各方向成分の和はそれぞれ独立に計算することができるため，直交格子法の幾何学的な特性を利用して，面の垂直方向と面積を容易に得ることで簡便な計算が実現できる．例として，x 方向には，nx = ±1 および ∆S =∆y∆z を用いることができる．他の方向に対しても同様である．ここで，必要となる面積分は，物体内部（ゴーストセル）と境界セルの境界面だけとなる．
次に，流体力の評価に必要な格子面での中点（Fig. 3(b)の Ix と Iy ) で応力の求め方を示す．ここでは 2 種類の方法を示す．まず 1 つ目の方法（MM1) に関して示す．本手法は，後述する 2 つ目の方法 (MM2) の技術のキーを示すため，あえて不適切な応力の評価方法を行っている．
MM1 における圧力の評価方法は，

pI,xpI,yI1

facet

(12)(13)

face,yIyvoxel)face)center
boundary)cellsolid)cell(b) Mesh-based method

Fig. 3: Polygon and mesh based methods.とした．ここで pI,x と pI,y はそれぞれ，Ix および Iy での圧力値を示している．また，せん断応力の近似に関しては，Ix および Iy での速度を物体速度 uI = 0 として，∂u∂x I,x=ui,j −uI∆x2,(14)

3c 2013 by JSFM
Copyright ⃝第 27 回数値流体力学シンポジウム
C02-3∂u∂y I,xui,j+1 −ui,j−1,2∆x=(15)

ここで, 計算に用いたステンシルを Fig. 4 に示す．
次に 2 つ目の適切な評価による方法 MM2 の応力評価方法を示す．圧力に関しては，

pI,xui−1, j+1

pI,yui, j+1Ix∂u∂x I,xui−1, j−1(16)(17)

また，せん断応力の近似に関しては，2 次精度中心差分を用いた．点 Ix における速度の x,y 方向微分を以下のように算出した：

ui, jui−1, jpi,j +pi−1,j,2pi,j +pi,j−1=,2=∂u∂y I,xui, j−1ui,j −ui−1,j,
∆x=
= 21

[ui,j+1 −ui,j−1u−u+ i−1,j+12∆xi−1,j−12∆x(18)], (19)

本計算では,x, y 方向微分を行う際に物体内部の値を必要とする．この値は, 数値解析を行う際にも必要となるものである．本研究では，上述の Eqs. 6, 7 を用いて これらの値の算出に面の垂線方向に伸ばした線上のイメージポイントから外挿した． (Fig. 4),

ﬂuid%cell
boundary%cell
solid%cellpS = p(−d0 ),uS = u(−d0 ),(20)

ここで，d0 は物体表面から固体セルの中心までの距離であり， d1 と d2 はそれぞれ 1.75∆x と 3.5∆x に設定した．

(a) Stencil for the computation of the derivatives

Stokes 流れの流体力評価の精度
上記で導入した流体力評価方法の精度を議論するため,Stokes 流れの厳密解を用いてその流体力を評価する.

3.

probe

3.1 球まわりの Stokes 流れ
球周りの Stokes 流れの厳密解は以下のように与えられる．

I1pI2

I

d1 d2uv

Sd0w

3µaxu,
2 ρ ∞ r3()()
1a23a2aax2= u∞ − u∞3 + 2 − u∞ 3 1 − 2 ,4rr
4rr()
23axya= − u∞ 3 1 − 2 ,
4rr()
3axz
a2= − u∞ 3 1 − 2 ,(21)4rr= p∞ −

ここで a は球の半径，ρ は密度，µ は粘性係数，u∞ は一様流速度，p∞ は参照圧力である．また，x，y および zはデカルト座標であり，r は球の中心にとった原点からの距離である．
流れの方向は x 方向を仮定しており，球の抵抗は Eq. 1で定義される力の x 方向成分となる．この定義式に厳密解である式 21 を代入すると，圧力抵抗と粘性抵抗がそれぞれ与えられる．

solid'cell
points'for'interpola-on

(b) Interpolation for the solid cells

2πaµu∞Fxp8µ8===, (22)2 /2) (πa2 )2 A(ρuρuDReρu∞∞∞2Fv
4πaµu∞
16µ16
CDf = 1 x2===, (23)2 /2) (πa2 )(ρuρuDReρuA∞∞∞2

CDp = 1Fig. 4: Stencils for evaluation of stresses and interpolation to the solid cell.

ここで，CDp と CDf はそれぞれ圧力抵抗係数と粘性抵抗係数であり，Re = ρu∞ D/µ はレイノルズ数, D = 2aは球の直径，A = πa2 は球の射影面積となっている．

4c 2013 by JSFM
Copyright ⃝第 27 回数値流体力学シンポジウム
C02-33.2 計算条件
本節では，厳密解を用いるため，特に計算領域を広くとる必要はないが，ゴーストセルの取り扱いがあるため，計算領域を L = 2D とし，中心に球を配置した Stokes 流れを厳密解から得た．この計算領域を 2N のセルに分割し，解析を行った．格子点は，N = 4, 8, 16, 32, 64 および 128 と設定した．それぞれのイメージ図を Fig. 5 に示す．PM に関しては，ポリゴンが必要になるが，球座標系の等間隔の 400 × 400 のポリゴンを用いた時に，今回調べた全ての計算格子で力の収束が見られたため，これを用いた．
3.3 結果Fig. 6 に PM および MM1，MM2 で評価した圧力抵抗及び粘性抵抗を示す．図中に，800 および 1600 が厳密解である．PM，MM2 に関してはほぼ同じ精度で流体力を評価できることがわかる．このことから，今回提案する方法は，ポリゴンに内挿する PM とほぼ同程度の精度で流体力評価ができるといえる．MM2 の計算コストは，PMに比べ非常に低いこと，また計算コードが簡単に書けるため，見通しの良い負荷の低い実装が可能になる．一方で，MM1 は圧力抵抗に関しては収束しているが，粘性抵抗は収束していないことがわかる．
ここでエラーの切り分けを考える．MM1,MM2 共に持っているエラーとして，1) 積分面が実際の物体と異なるためのエラーおよび 2) 応力の評価におけるエラーがある．
まず後者のエラーを打ち消すため，応力の格子面の位置での厳密解を与えて解析したものを Fig. 7 に示す．前
述の結果と比較すると，厳密解を与えた MM-Analiticalは MM2 とほぼ同じ結果となっており，厳密解に収束していく．このことから 1) の積分面が実際の物体と異なるためのエラーは問題なく，2) 応力の評価におけるエラーを減らすことが重要といえる．すなわち，MM1 では，格子面上で速度 0 を置き間違った応力を計算で使ってしまうために流体力が正しく評価できないと言える．
1) に関しては，式からわかるように表面の面積ベクトルの積分値が正しい値に収束すれば良いと考えられる．
Fig. 8 に MM 法で用いた表面の面積積分値と面積ベクトルの近似値が格子解像度に対してどのように変化するかを示した．MM 法では，面積の積分値は直交格子を使っているため，厳密な値には収束しないが，面積ベクトルの成分は，厳密値に収束している．すなわち式で必要な面積ベクトルの各成分は，直交格子の格子面を用いても正しい値に収束させることができる．これらの結果は，流体解析で得られた数値結果にある誤差を取り除いた上で流体力の評価方法のみを議論することで初めて明らかに
5700

600500PMMM1MM2400300
0.001

0.010.1

10.1

1∆x(a) Pressure drag2000

Pressure drag coeffients

Fig. 5: Spheres expressed by the meshes with diﬀerentresolutions. Here, the sphere expressed by the mesh of∆x = a/4, x = a/8, a/16, a/32, a/64, a/128 are shown
from left to right.

Pressure drag coeffients

8001800

1600

1400
PMMM1
MM2

1200

1000
0.001

0.01

∆x(b) Friction dragFig. 6: Estimated drags for Stokes ﬂow as the functionof mesh spacings. Here, 800 and 1600 are analyticalestimation of pressure and friction dragsc 2013 by JSFM
Copyright ⃝x-component of surface area vectors

第 27 回数値流体力学シンポジウム
C02-3

Pressure drag coeffients

800

700

0.3
0.28
0.26
0.24
0.22
0.2
0.18
0.16
MMexact

0.14
0.12
0.1
0.001

0.01

500

1(a) x component of surface area vector

MM1
MM with analytical shear stress

400
0.7
300
0.001

0.65
0.01

0.1

1

∆xSurface area

0.6(a) Pressure drag

0.55
0.5
0.45
0.4

2000

0.35

Pressure drag coeffients

0.1

∆x

600

1800

0.3
0.001

0.01

0.1

1

∆x

1600

1400

1200

MMexact(b) surface area

MM1
MM with analytical shear stressFig. 8: Evaluation of surface area and surface area vector for MM.
1000
0.001

0.01

0.1

1

∆x
出来た．

Re=100 の球流れ
ここでは，実際に Re=100 の球周りの流れを別途開発した直交格子ソルバー(10) で，解析した結果を用いて前節と同様に，流体力の評価を行った．本ソルバーは,SMAC 法を基礎とし, 佐藤ら(2, 3) と同様の Immersed boundary 法の実装となっている. 本評価では, 数値解析そのもののエラーも含まれていることに注意されたい. また, 格子点数の増加を防ぐため, 階層型直交格子(11) を採用している．球付近の格子幅,∆x/2a = 0.05859375 ∆x/2a = 0.029296875∆x/2a = 0.0146484375 の 3 ケースを用いて，評価を行った．このレイノルズ数では，抵抗係数 1.08 程度になることが知られている．Fig. 9 に示すように，MM2 は PMと同様の精度で流体力を予測できることがわかる．

4.

(b) Friction dragFig. 7: Error decomposition of drags for Stokes ﬂow
using the analytical shear stress. Here, 800 and 1600
are analytical estimation of pressure and friction drags

結論
本論文では,Immersed boundary 法における格子をベースにした簡便な流体力の見積もり方法を提案した．Stokes

5.

6c 2013 by JSFM
Copyright ⃝第 27 回数値流体力学シンポジウム
C02-3(6) Yamazaki, W., Matsushima, K., and Nakahashi, K., “Aerodynamic design optimization using
the drag-decomposition method,” AIAA journal ,Vol. 46, No. 5, 2008, pp. 1096–1106.

1.15
1.14

Drag coeffients

1.13(7) Lai, M.-C. and Peskin, C. S., “An immersed boundary method with formal second-order accuracy andreduced numerical viscosity,” Journal of Computational Physics, Vol. 160, No. 2, 2000, pp. 705–719.

1.12
1.11
1.1

PMMM1
MM2

1.09

(8) Balaras, E., “Modeling complex boundaries using
an external force ﬁeld on ﬁxed Cartesian grids
in large-eddy simulations,” Computers & Fluids,Vol. 33, No. 3, 2004, pp. 375–404.

1.08
1.07
1.06
1.05
0.01

(9) Shen, L., Chan, E.-S., and Lin, P., “Calculationof hydrodynamic forces acting on a submerged
moving object using immersed boundary method,”Computers & Fluids, Vol. 38, No. 3, 2009, pp. 691–702.

0.1∆xFig. 9: Estimated drags for Reynolds number 100 as
the function of mesh spacings.

流れの厳密解を用いることで，流体解析に含まれるエラーを排除して，流体力の見積もり方法を評価した．本論文で提案する方法は，ポリゴンをベースにした方法と同程度の精度を持ちながら，その計算コストは非常に低く，簡便な実装が可能である．
この手法を評価する上で，流体力評価のエラーとして，積分面が格子面となっていること，および応力の評価方法の 2 つがあると考えそれぞれを検討した．応力の評価方法を適切に行うことで，積分面が格子面となっていても，本手法は十分な精度が得られることがわかった．これは，格子面はどんなに細い格子を用いても，その積分面積は物体の表面積分値に一致しないが，一方で，流体力評価に現れる物体の表面ベクトルの積分値はよく一致するためである．また応力の適切な評価は，物体内部点を用いた差分公式を使うことで実現できることがわかった．

(10) Onishi, J., Ono, K., and Suzuki, S., “Developmentof a CFD software for large-scale computation: Anapproach to grid generation for arbitrary complex
geometries using hierarchical blocks,,” to be published in Proceedings of the 12th International Symposium on Fluid Control, Measurement and Visualization(FLUCOME2013), 2013 , 2013.
(11) Nakahashi, K., “High-density mesh ﬂow computations with pre-/post-data compressions,” AIAAPaper 2005-4876 , 2005.

謝辞
第一著者は科研費 (25709009,24656522) の支援を受けた．また本研究は，HPCI 戦略プログラム 分野４の支援を受けた．本稿に対して，JAXA 青野光博士に貴重なコメントを頂いた．ここに記して謝意を表する．
参考文献(1) Mittal, R. and Iaccarino, G., “IMMERSEDBOUNDARY METHODS,” Annual Review ofFluid Mechanics, Vol. 37, 2005, pp. 239–261.
(2) 佐藤範和, 竹内伸太郎, 梶島岳夫, 稲垣昌英, and 堀之内成明, “直交格子を用いた対流熱伝達計算における物体境界近傍の直接離散化法,” 日本機械学会論文集 B 編, Vol. 79, No. 803, 2013, pp. 1219–1231.
(3) 佐藤範和, 梶島岳夫, 竹内伸太郎, 稲垣昌英, and 堀之内成明, “直交格子法における物体境界近傍の直接離散化法 (速度場と圧力場の整合性を考慮した高精度化),” 日本機械学会論文集 B 編, Vol. 79, No. 800,
2013, pp. 605–621.
(4) 高橋悠一 and 今村太郎, “直交格子法を用いた粘性計算における力計算と物体壁面境界の取扱いについて,” 第 26 回数値流体力学シンポジウム講演集, 2012.
(5) Takahashi, S., Nonomura, T., and Fukuda, K.,“High-resolution Immersed Boundary Method forCompressible Turbulent Flows with Shocks: Application to Two-dimensional Flows around Cylinders,” submitted., 2013.
7c 2013 by JSFM
Copyright ⃝