オイラー型統一解法による固体-流体連成解析手法
中尾賢司 1) ，岡澤重信 2) ，西口浩司 3) ，杉山和靖 4) ，野田茂穂 5) ，高木周 6) ，岡田裕 7) ，松澤照男 8) ，熊畑清 9) ，北脇知己 10) ，天野晃 11)
1)

広島大学 大学院工学研究科，nakao-kenji@hiroshima-u.ac.jp
広島大学 大学院工学研究科，okazawa@hiroshima-u.ac.jp
3)
広島大学 大学院工学研究科，nishiguchi@hiroshima-u.ac.jp
4)
東京大学 大学院工学系研究科，sugiyama@fel.t.u-tokyo.ac.jp
5)
理化学研究所 情報基盤センター，shigeho@riken.jp
6)
理化学研究所 臓器全身スケール研究開発チーム，takagish@riken.jp
7)
鹿児島大学 大学院理工学研究科，okada@mech.kagoshima-u.ac.jp8)
北陸先端科学技術大学院大学 情報科学センター，matuzawa@jaist.ac.jp
9)
北陸先端科学技術大学院大学 情報科学研究科，k kuma@jaist.ac.jp
10)
岡山大学 大学院保健学研究科，kitawaki@md.okayama-u.ac.jp
11)京都大学 大学院情報学研究科，amano@i.kyoto-u.ac.jp
2)

生体力学シミュレーションのための固体-流体連成解析手法を提案する．本手法ではボクセルデー
タと親和性の高いオイラー型解法をベースとしてスキームを構築する．固体と流体の双方を非圧縮性の物質と仮定し，それぞれの密度関数 (要素内体積率) を用いて両者の支配方程式を統一化する．そして固体-流体連成問題に適用することで，本手法の妥当性を検証する．

1.

はじめに製造分野では通常，解析対象のモデルを CAD などを用いて作成し，計算対象のモデルとして利用する．一方，生体力学シミュレーションの分野では，CT や MRI等の画像診断装置を用いて解析対象の形状を計測し，そこから得られるデータを計算対象のモデルとして利用する．このデータはボクセルと呼ばれる形式であり，計測した値の空間分布を示している．一般的に，このデータは従来の製造業で用いられてきた力学シミュレーションとは整合性が悪いため，ボクセルデータを直接シミュレーションの中に取り込み，効率的に解析を行うような手法の開発が求められる．
一般的に，ボクセルデータに適した計算方法としてはオイラー型解法が挙げられる．オイラー型解法は空間に固定された座標から物体の変形を評価する方法であり，通常は流体の解析に用いられることが多い．本
研究ではオイラー型解法をベースにスキームの構築を目指す．
生体は本来柔らかく，血流の変化などの影響で容易に変化するため，その再現には固体-流体連成解析が必要とされる．従来の連成解析では，図-1(a) のようにオイラー的に記述された流体に対して，固体はラグランジュ法により表現して解析を行う方法が一般的であった．しかし，この方法ではボクセルデータと親和性が悪く，また固体側のメッシュの破綻により，生体の大変形を表現できないという問題が発生する．

以上のような背景を踏まえ，本研究では図-1(b) のよ
うに，固体に対してもオイラー的に記述し，固体と流体をオイラー表記により統一的に表したボクセル解析のための固体-流体連成解析手法の構築を目指す．固体と流体の運動記述を支配方程式の段階で統合し，統一的解法として解析を行う方法を提案する．

(a) Conventional approach(b) Present approach

図–1 Conventional and present solid-ﬂuid interaction approaches2.

固体-流体統一的解法

本研究では，固体を非圧縮性弾性体，そして流体を非圧縮性流体と仮定する．今，図-2 のような固体と流体の双方が存在する検査体積を考える．なお，この検査体積の体積を V ，表面を S で表している．
図-2 の固体領域では以下の支配方程式が成り立つ．

ρsDus= ∇ · σs + ρs bDt

(1)

S

2.1 構成方程式
非圧縮性を仮定した場合，弾性体とニュートン流体の応力は一般的に以下のように偏差成分と圧力成分に分けた式で表すことができる．

FluidVSolid

図–2 Solid-Fluid model∇·us = 0(2)ここで，添え字の s は固体の物理量であることを示している．また，流体領域では以下の支配方程式が成り立つ．
Duf= ∇ · σf + ρf b(3)ρfDt∇·uf = 0(4)ここで，添え字の f は流体の物理量であることを示す．
それぞれの領域で，それぞれの支配方程式が成り立っているため，検査領域全体では以下の支配方程式が成り立つ．
Duf
Dus+ If ρfIs ρs
Dt
Dt= Is ∇ · σs + If ∇ · σf + (Is ρs + If ρf )b(5)Is ∇ · us + If ∇ · uf = 0(6)3.

(ϕs ρs + ϕf ρf )D(ϕs us + ϕf uf )Dt

UWFluid(9)∇ · (ϕs us + ϕf uf ) = 0(10)ここで，上付きのバーは領域全体として共通の物理量であることを表している．この時，支配方程式 (5) の体積平均は，数学的には直接式 (9) の形にはならない．
解析上取り扱いが容易な式 (9) の形にするために，省略している項がある．この項は平均化，混合化により派生する項であり，省略する根拠としては空間解像度を高くすると，これらの影響が小さくなることが挙げられる．この項については紙面の都合上，詳細を割愛する．

0.50.5Solid0.50.5-UW

図-3 Solid-Fluid shear test

図-4 Velocity proﬁle

図-4 は解析結果であり，周期的に変化する速度分布を表している．流体層と固体層で速度の位相が変化する様子が捉えられている．この結果は本研究の統一化
支配方程式から導いた理論解ともよく一致している．

4.

= ∇ · (ϕs σ s + ϕf σ f ) + (ϕs ρs + ϕf ρf )b(12)

本手法での固体-流体連成解析の妥当性を検討するために，図-3 に示すような，流体層に挟まれた固体をせん断運動させる問題を考える．このとき，無次元量で
密度 ρs = ρf = 1，粘性係数 µ = 1，せん断弾性係数 G = 5 としている．また上下境界に与える速度は
UW = ÛW sinωt としている．

(7)以上の関係を用いて，支配方程式 (5),(6) を体積平均し，以下のような形の式で表す．

σf = σf′ − pf I

解析結果

Fluid

ここで，⟨ ⟩ は体積平均操作を意味する．また，ϕs ,ϕfは固体，流体それぞれの体積率を表し，以下の関係を満たす．
ϕs + ϕf = 1(8)(11)

ここで，偏差成分 σs′ と σf′ は，固体はフック弾性体の構成式，流体はニュートン流体の構成式を用いて計算することができる．圧力に関しては固体，流体双方とも非圧縮拘束条件よりポアソン方程式を導入して求める(SMAC 法を適用)．この時，固体と流体が混ざり合った要素では，平均された圧力 pmix = ϕs ps + ϕf pf が直接求まる．

ここで Is ,If は Indicator 関数であり，例えば Is は固体領域では 1，その他の領域では 0 の値を取る．ここで，支配方程式 (5),(6) を体積平均することで，平均化方程式を導く．Indicator 関数 Is ,If は体積平均を行うと以下のようになる．

⟨Is ⟩ = ϕs , ⟨If ⟩ = ϕfσs = σs′ − ps I

まとめ

本研究では，ボクセルデータと親和性の高いオイラー型解法を用い，固体-流体連成解析手法の構築を目指した．固体と流体それぞれの体積率により支配方程式を統一し，統一的解法の提案を行った．連成計算を行った結果，妥当な解が得られ本手法の有用性を示すことができた．
参考文献1) 梶島岳夫 : 乱流の数値シミュレーション，第 1 章-第 3章，養賢堂，1999．
2) Benson D.J.：Computational methods in Lagrangianand Eulerian hydrocodes, Comput．Methods Appl．
Mech．Engrg., 99, 235-394，1992.