第 24 回数値流体力学シンポジウム
B5-4

血管異方性を考慮した超弾性体管内流のオイラー型流体構造連成シミュレーション
Full Eulerian fluid-structure coupling simulations of fluid flow in anisotropic hyperelastic tube⃝

長野 直大，東大院，〒 113-8656 東京都文京区本郷 7-3-1，E-mail：nagano@fel.t.u-tokyo.ac.jp
杉山 和靖，東大工，〒 113-8656 東京都文京区本郷 7-3-1，E-mail：sugiyama@fel.t.u-tokyo.ac.jp
竹内 伸太郎，阪大工，〒 565-0871 大阪府吹田市山田丘 2-1，E-mail：shintaro@mech.eng.osaka-u.ac.jp
伊井 仁志，東大工，〒 113-8656 東京都文京区本郷 7-3-1，E-mail：sii@fel.t.u-tokyo.ac.jp
高木 周，東大工（理研），〒 113-8656 東京都文京区本郷 7-3-1，E-mail：takagi@mech.t.u-tokyo.ac.jp
松本 洋一郎，東大工，〒 113-8656 東京都文京区本郷 7-3-1，E-mail：ymats@mech.t.u-tokyo.ac.jp

Naohiro NAGANO, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-8656
Kazuyasu SUGIYAMA, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-8656
Shintaro TAKEUCHI, Osaka University, 2-1, Yamadaoka, Suita-shi, Osaka, 565-0871
Satoshi II, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-8656
Shu TAKAGI, The University of Tokyo(Riken), 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-8656
Yoichiro MATSUMOTO, The University of Tokyo, 7-3-1, Hongo, Bunkyo-ku, Tokyo, 113-8656
An interaction problem between a fluid and elastic walls is solved by a full Eulerian fluid-structure coupling
method. The method employs a uniform grid system for both fluid and solid and it does not require any meshgenerarion or reconstruction, aming for facilitaitng the practical bio-mechanical fluid-structure analysis. The
availability and applicability of the developoed method to systems involving complex geometries driven by a pressure gradient are shown through a comparison between the obtained numerical results and theoretical prediction,
and a grid convergence test. Further, for considering a more realistic situation, an anisotropic blood vessel modelis implemented and convergence performance of the method is shown.

はじめに
近年，医療応用の観点から，粘性流体中における固体の変形・運動を扱った流体構造連成シミュレーションが注目されている．特に，CT や MRI といった画像診断機
器の技術的向上に伴い，人体内部の情報をより正確かつ，
簡単にボクセルデータとして取得可能になってきたことから，患者個人のボクセルデータに対して流体構造連成
シミュレーションを行うことで，各人の治療効果の予測や治療方針の決定に役立てられると期待されている．さらに，患者の全身ボクセルデータを用いた流体構造連成
シミュレーションを行うことで，患者個人にテーラーメ
イド型医療を提供することが期待されている．
生体の流体構造連成シミュレーションでは，生体内の複雑な形状に適応するため，固体の変形をラグランジュ
的に記述する有限要素法 (FEM) がもっぱら用いられてきた (1)−(3) ．FEM による解析では，CT や MRI から得られたボクセルデータを元にメッシュを生成し，時間ス
テップごとに，メッシュの再構成が必要となる．患者個
人のボクセルデータを用いた FEM による流体構造連成
シミュレーションとしては，脳動脈瘤付近の血流の解析(4),(5)などがある．
また，Matsunaga や Yokoi らは血管壁を剛体壁とし，脳動脈瘤付近の流れの解析を診断画像から得られたデー
タを直接扱える手法で行ってきた (6),(7) ．
一般に，血管壁はいくつかの異なる力学的性質を示す層からなることが知られており，血管内膜の変性や崩壊が脳動脈瘤の発生原因となる可能性，血流によってさらに変性が進行し，脳動脈瘤となる可能性が指摘されている (8),(9) ．そのため，血管壁の変形を考慮した血流シミュ
レーションが脳動脈瘤といった血管が一定の形状を保たない場合において重要となる．
これらの理由から，本研究では，血管の変形を考慮し，かつ診断画像から得られるボクセルデータを直接扱える
血流シミュレーション手法の開発を目的としている．これまでに，杉山らが開発した，オイラー型定式化を用いた

1.

1差分法による流体構造連成シミュレーション手法 (18),(19)を拡張し，血管壁を等方性の超弾性体として模擬し，本計算手法の妥当性の検証などを行ってきた (21) が，一般に，血管壁はその力学的性質が異方的であることが知られており，様々な血管壁モデルが提案されている (10)−(15) ．
本論文では，これまで開発してきた計算手法の概要とその検証，さらに血管異方性の導入，その検証について述べる．

2.

数値計算モデル

2.1

支配方程式

本研究では，多くの生体に対する流体構造連成シミュ
レーション (16),(17) と同様，非圧縮性の流体・固体を仮定
し，流体・固体の密度は同じとしている．流体・固体は以下の質量保存式：∇ · v = 0,(1)

運動量保存式：(ρ∂v+ v · ∇v∂t)

= −∇p + ∇ · σ ′ −

∆P
ex
L

(2)

に従う．ここで，v(= (vx , vy )) は速度ベクトルを，ρ は密度，p は駆動圧力からの圧力の変動分，σ ′ は偏差コー
シー応力テンソル，−∆P/L は駆動圧力，ex は x 方向への単位ベクトルである．また，これ以降 ′ はテンソルの偏差成分を表す．

c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
B5-4

2.2 構成方程式
流体・固体を区別するために，VOF(Volume-of-Fluid)
法のアイディアに基づいた関数 ϕs を導入する．ϕs は
1 (solid)(3)ϕs =0 (ﬂuid)行う．なお，Holzapfel らは，血管壁全体の力学的性質として Intima がほとんど効果を持たないとして，血管壁をMedia，Adventitia の 2 層からなる固体として 3 次元モデルの構築を行っているが，本論文では，まずは 1 層の
2 次元モデルとして導入を行う．
オイラー場で固体の応力を表現するために，初期状態の繊維の単位ベクトル a0 ，時間 t における繊維のベクト
ル a(= F · a0 )，さらに a を用いて，となる固体の体積率を示すスカラー関数である．数値計算
上では，流体・固体を単一の連続体と捉え，偏差 Cauchy
応力 σ ′ を ϕs を用いて各相の応力を体積平均することで
定義する：
(4)
σ ′ = (1 − ϕs ) σf′ + ϕs σs′ .

Ha = a ⊗ a = F · (a0 ⊗ a0 ) · F Tを定義する．Ha は現配置によって定義されるテンソルのため，Ha を用いて固体の応力を以下のようにオイラー場で表現できる：

ここで，添え字 f は流体を，s は固体を示す．流体のモデ
ルとしてニュートン流体，固体のモデルとして neo-Hooke
粘弾性体を用いる．流体・固体の応力はそれぞれオイラー場において，以下のように表せる：

σfσs= 2µD,= 2µD + GB.

σs = 2µD + GB()
2
+ 2k1 exp k2 (Ia − 1) (Ia − 1) Ha . (12)

(5)
(6)

ここで，k1 ,k2 は定数，Ia = trHa である．
Ha は B 同様，Oldroyd 速度が恒等的に零となる物理量であり，以下の輸送式に従う：

ここで，µ は粘性係数，D(= 12 (∇v + ∇v T )) は速度勾配テンソル，G は横弾性係数，B(= F · F T ) は左 CauchyGreen 変形テンソル (F (= ∂x/∂X) は変形勾配テンソル)である．
また，物理量 ϕs ，B は以下の輸送式に従う：∂ϕs+ v · ∇ϕs = 0,∂t∂Ha
+ v · ∇Ha = L · Ha + Ha · LT .
∂tH̃a = ϕ1/2
s Ha

(7)

(14)

を定義し，実装上では，偏差 Cauchy 応力 σ ′ は′
σ ′ = 2µD ′ + ϕ1/2
s GB̃
(
)
2
+ 2ϕ1/2(Ia − 1) H̃as k1 exp k2 (Ia − 1)(15)としている．これ以降，偏差 Cauchy 応力が式 (15) に表される計算モデルを異方性モデルと呼ぶ．

(9)を定義し，実装上では，偏差 Cauchy 応力 σ ′ は′
σ ′ = 2µD ′ + ϕ1/2
s GB̃

(13)

また，B 同様，計算安定化のため，∂B+ v · ∇B = L · B + B · LT .
(8)∂t
ここで，L(= ∇v) は速度勾配テンソルである．式 (8) は，
物理量 B が恒等的に Oldroyd 速度が零となることを示す．
また，杉山らにならい (18),(19) ，計算安定化のため，
B̃ = ϕ1/2
s B

(11)

(10)

としている．これ以降，偏差 Cauchy 応力が式 (10) に表される計算モデルを等方性モデルと呼ぶ．

2.3 血管異方性の導入
本研究では，血管のオイラー型流体構造連成シミュレー
ション手法の基礎として，血管壁を等方的な超弾性体としてモデル化して手法を開発してきた (21) ．しかしながら，一般に，血管壁は異なる力学的性質を持つ Intima，Media，Adventitia の 3 層から構成されること，力学的性質として異方性を示すことなどが知られており，本論
文では血管異方性を考慮した材料モデルの導入を行う．
血管の材料モデルとしては，Fung らが提案したモデル(11)や，それに基づいたモデル (13) などが良く知られている．Fung らは指数関数を用いて歪みエネルギー関数を定義し，in vitro での実験値と良く合うモデルの構築を行った．Fung らが構築した材料モデルは一様な異方性材料からなるのに対し，Holzapfel らは，血管壁を構成する
層ごとに，等方的な超弾性体と，異方性を持つらせん状の繊維によって構成される材料モデルを考え，定式化を行った (15) ．なお，Fung らのモデルと同様，異方性の成分については指数関数を用いて歪みエネルギー関数を定義している．
本論文では，2.2 節に記した今まで開発してきた計算モ
デルを拡張し，Holzapfel らの異方性血管モデルの導入を
2

2.4 計算条件
図 1 に示したような，弾性体壁を含む波状流路内の非定常ストークス流れの解析を行う．そのため，式 (2) の対流項は無視しているが，式 (7),(8),(13) の対流項はそれぞれの物理量の時間発展を追うため，考慮している．
計算領域は上下を剛体壁に囲まれており，弾性体と剛体壁の境界において no-slip 条件を課している．そのため，その境界上において，式 (8) の各成分は，以下のよ
うに展開できる：B̃xx2= 1 + B̃xy
,

B̃yy

=

ϕ0s = 1,∂ B̃xy
∂t

=

∂vx
.
∂y

(16)
(17)
(18)

ここで，添え字 0 は初期値を示す．また，式 (16) の導出に非圧縮性材料の条件である，detB = 1 を用いた．
同様に式 (13) の各成分は，以下のように展開できる：∂ H̃axx
∂t

=

2

∂vx
H̃axy ,
∂y
0

H̃ayy

= H̃a yy ,∂ H̃axy∂t= H̃a yy

0

∂vx
.
∂y

(19)
(20)
(21)

c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
B5-4

また，x 方向に波長 L の周期境界条件を課している．
図 1 の変数は本研究では H0 = h0 = 1，k = π/4，
L = 8，δ = 0.3 としている．
初期条件として速度場と圧力場は静的平衡状態を仮定
し，ϕs を各格子ごとに計算して与える．また初期状態において，変形勾配テンソル F は，

このとき基底関数の直交性から，フーリエ係数は以下のように表せる：
(
)
(
)
∫
Hcn
cos 2πnx
2 L
L=dx H(x)
.
(27)
L 0
Hsn
sin 2πnx
L

F =I

ここで，固体体積率関数 ϕs がヘビサイド関数状に変化することを利用し，フーリエ係数は，以下のように表せる：
(
)
(
)
∫
∫ L
Hcn
cos 2πnx
2 H0 +h0
∂ϕs
L
≃−
.
dy
dx
y
L 0
∂y
Hsn
sin 2πnx0L(28)
離散化された上式は次式のようになる：(22)となるため，左 Cauchy-Green 変形テンソル B の初期値は
B = I,(23)

Ha の初期値は
Ha = I · (a0 ⊗ a0 ) · I

(24)

(

として与えられる．また，初期状態の繊維の単位ベクト
ル a0 を



1

√




 1+δ2 k2 cos2 ( 2πxL )a0 = 
2πx


−δk
cos
( L )  (0 ≤ y ≤ 2)

√


2 2
2 2πx
 1+δ k cos ( L ) 
(25)
1

√



 1+δ2 k2 cos2 ( 2πx


L )

a0 = 2πx (2 ≤ y ≤ 4)

δk
cos
( L )

√


1+δ 2 k2 cos2 ( 2πx
L )

)

2 ∑∑
(ϕs (i, j + 1) − ϕs (i, j))
L i j
(
)
cos 2πnx(i)Ly(j + 1/2)(∆x). (29)sin 2πnx(i)L≃−

得られたフーリエ係数を元に以下のようにスペクトルを定義する：√
2 + H2 .
XnNx = Hcn(30)sn

G = 2 の場合について，最高解像度 (Nx × Ny = 1024 ×
512) の計算結果から得られるスペクトルと，他の解像度から得られるスペクトルとの相対誤差を図 4 に示す．図4 から，低次のモードに関して，相対誤差は格子解像度にほぼ比例して減少していることが分かり，弾性体壁が剛体壁に固着し，歪みが蓄積されていく本計算条件下においても，格子解像度に対してほぼ一次の誤差精度を有していることが分かる．

と与える．これは，図 2 に示すように界面に沿った繊維を仮定している．
また，本論文を通じて，t = 0 以降，一定の圧力勾配−∆P/L (= 1) を系全体の駆動力として与える．

2.5 数値計算法
全ての方程式は，矩形固定スタッガード格子上で差分法により解かれる．空間微分の離散化については式 (7)と (8)，(13) の対流項には五次精度の WENO 法を，それ以外の項には二次精度の中心差分法を用いる．時間進行法については，偏差 Cauchy 応力テンソルには二次精度の Crank-Nicolson 法を，それ以外の項については二次精
度の Adams-Bashforth 法を用いる．速度場のソレノイダ
ル条件を満たしながら，運動量保存式の時間進行を行う
ために，非圧縮性流体の標準的解法である SMAC 法を用いる．また，x 方向に高速フーリエ変換を用いて圧力のPoisson 方程式を解く．
3.

Hcn
Hsn

3.2 妥当性の検証
図 1 のように流体・弾性体壁の界面が波状となっている問題を対象として，本計算条件下における本計算手法の妥当性の検証を行う．本計算条件下においては，系全体は一様な圧力勾配で駆動されており，一定時間経過後，系全体は定常状態に達する．そこで，定常状態時の界面位置の理論解と得られた計算結果を比較することで，妥当性を検証する．
本研究では大変形を伴う超弾性体を扱っており，弾性体の変形に伴う幾何学的非線形影響が，一般には無視できないため，厳密解を求めることは困難である．本研究では，流体相でレイノルズ数が十分に小さく，長波近似が成り立ち，また弾性体相で線形弾性体とみなせる程度に変形が微小となる条件において，定常状態における弾性体の変位量を潤滑理論と組み合わせた線形解析により求めた．
図 5 に G = 10 の場合について，流れ場が十分に発達し，ほぼ定常状態に達したと考えられる，t = 30 の超弾性体壁の変形形状を固体体積率関数 ϕs のコンターとして示す．また，比較のため超弾性体壁の初期形状を点線で示す．なお，格子点数は x，y 方向それぞれについて
Nx = 256，Ny = 128 であり，時間解像度は ∆t = 0.001とする．
先に得られた理論解 (21) と，計算結果から得られた界面位置を比較することで，妥当性の検証を行う．図 6 に比較した結果を示す．得られた理論解と，計算結果から得られた流体と弾性体壁の界面位置は，図 6 に示される
ように良く一致する．これは理論解を導く上で用いた仮定が G = 10 の弾性体によく合っていたことを示すと同時に，本計算手法が本計算条件下においても妥当である
ことを示している．以上から本計算手法が圧力勾配により駆動される複雑形状を有する流れ場に対しても十分な精度で適用可能なことが分かる．

等方性モデルの検証

3.1 収束性の検証
杉山らは超弾性体粒子を含む流れ場の解析を行い，固体の材料モデルとして neo-Hooke 体や線形 Mooney-Rivlin
体，St.Venant-Kirchhoﬀ 体を用いて，本計算手法が格子解像度に対して一次の誤差精度を有していることをすでに示している．しかしながら，計算条件が異なるため，本計算条件下における格子解像度に対する誤差精度を確認する．
固体体積率関数 ϕs は界面において不連続な値の分布を持つため，陽に界面位置を定義することができない．そのため，流れ場が十分に発達したと考えられる t = 30 の計算結果を用いて，以下のような操作を行い，界面位置を定義する．なお，t = 30 における流れ場の様子を図 3に示す．
まず，底面からの弾性体壁の高さを H (x) とし，以下のようにフーリエ展開した形で表す：
)
∞ (
∑
2πnx
2πnx
Hcn cos+ Hsn sin.
H(x) = H0 +
λ
λ
n=1
(26)
3

c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
B5-4

異方性モデルの検証
ここから，2.3 節において述べた，弾性体の応力として血管異方性を導入したモデルを用いる場合について，計算手法の検証を行う．
まず，3.1 節と同様に界面位置をフーリエ解析により求め，スペクトルを定義，同様に最高解像度の計算結果との相対誤差を図 7 に示した．なお，計算条件は G = 2，
k1 = 1，k2 = 1 とし，十分流れ場が発達したと考えられ
る t = 30 の計算結果を用いた．このときの流れ場の様子を図 8 に示す．図 7 に示されるように，低次のモードに関して，相対誤差が格子解像度にほぼ比例して減少して
おり，3.1 節と同様，応力が式 (15) と表現される血管異方性を導入した際にも，格子解像度に対してほぼ一次の誤差精度を有していることが分かる．
また，式 (15) において係数 k1 を小さくすることで，計算結果は等方性モデルの計算結果に近づいていくと考えられる．そこで，十分流れ場が発達したと考えられる
t = 30 の計算結果を用いて，k1 の値を変化させ，界面位置を比較した図 9 に示す．図 9 から分かるように，k1 の値を減少させることで，等方性モデルの結果に近づいていくことが分かる．
以上の結果から，異方性モデルに関しても，本計算手法の収束性を確認できる．
なお，図 9 の x = 8 付近に見られる界面の振動の要因
として，主にせん断応力による弾性体の振動や，VOF 関数の界面表現精度法などが挙げられる．これらの問題は，著者らが別途取り組んできているシャープな界面追跡法(20)と陰解法 (20) との組み合わせにより高精度に解析できると考えられる．

4.

まとめ
既存のオイラー型流体構造連成手法を拡張し，血管壁を模擬した超弾性体管内流の数値解析を行った．圧力勾配により駆動される，超弾性体壁に囲まれたストークス
流れ場の流体構造連成問題を取り扱った．血管壁を等方的な超弾性体とした計算結果からは，定常状態の界面位置から収束性を確認し，定常解と理論解の良い一致から妥当性を確認した．さらに，Holzapfel らのモデルを基に，血管異方性の導入を行い，収束性の検証を行った．
今後は，さらに異方性モデルの検証を進め，三次元モデルの開発を行っていく．

5.

参考文献
(1) Taylor, C.A., Hughes, T.J.R. and Zarins, C.K.,”Finite element modeling of blood ﬂow in arteries”, Comput. Methods Appl. Mech. Engrg.
Vol.158(1998), pp.155-196.

(6) Matsunaga, N., Liu, H. and Himeno, R., ”An
immersed-based computational ﬂuid dynamics
method for haemodynamic simulation”, JSME Int.
J. Ser. C Vol.45(2002), pp.989-996.
(7) Yokoi, K., Xiao, F., Lui, H. and Fukasaku, K.,”Three-dimensional numerical simulation of ﬂows
with complex geometries in a regular Cartesian
grid and its application to blood ﬂow in cerebral artery with multiple aneurysms”, J. Comput.
Phys. Vol.202(2005), pp.1-19.
(8) Stehbens, W.E. ”Etiology of intracranial berry
aneurysms”, J. Neurosurg Vol.70(1989), pp.823831.
(9) Mizutani, T., Kojima, H., Asamoto, S. and Miki,
Y. ”Pathological mechanism and three dimensional
structure of cerebral dissecting aneurysm”, J. Neurosurg Vol.94(2001), pp.712-717.
(10) Vaishnav, R.N., Young, J.T. and Patel, D.J., ”Distribution of stresses and of strain-energy density
through the wall thickness in a canine aortic segment”, Circ. Res. Vol.32(1973), pp.577-583.
(11) Fung, Y.C., Fronek, K. and Patitucci, P., ”Pseudoelasticity of arteries and the choice of its
mathematical expression”, Am.
J. Physiol.
Vol.237(1979), pp.620-631.
(12) Takimizawa, K. and Hayashi, K., ”Strain energy
density function and uniform strain hypothesis forarterial mechanics”, J. Biomech. Vol.20(1987),
pp.7-17.
(13) Humphrey, J.D., ”Mechanics of arterial wall: Review and directions”, Critical Reviews in Biomed.
Engr. Vol.23(1995), pp.1-162.
(14) Delﬁno, A., Stergiopulos, N., Moore, J.E. and Meister, J.-J., ”Residual strain eﬀects on the stress ﬁeld
in a thick wall ﬁnite element model of the human
carotid bifurcation”, J. Biomech. Vol.30(1997),
pp.777-786.
(15) Holzapfel, G.A., Gasser, T.C. and Ogden, R.W.,”A new constitutive framework for arterial wall mechanics and a comparative study of material models”, J. Elasticity Vol.61(2000), pp.1-48.

(2) Torii, R., Oshima, M., Kobayashi, T. and Takagi,
K., ”Numerical simulation system for blood ﬂow inthe cerebral artery using CT imaging data”, JSMEInt. J. Ser. C Vol.44(2001), pp.982-989.

(16) Liu, W.K., Liu, Y., Farrell, D., Zhang, L., Wang,
X.S., Fukui, Y., Patankar, N., Zhang, Y., Bajaj C.,
Lee J., Hong J., Chen X. and Hsu H., ”Immersed
ﬁnite element method and its applications to biological systems”, Comput. Methods Appl. Mech.
Engrg. Vol.195(2006), pp.1722-1749.

(3) Zhang, Q. and Hisada, T., ”Analysis of ﬂuidstructure interaction problems with structural
buckling and large domain changes by ALE ﬁnite
element method”, Comput. Methods Appl. Mech.
Engrg. Vol.190-48(2001), pp.6341-6357.

(17) Huang, W.-X. and Sung, H.J., ”An immersedboundary method for ﬂuid-ﬂexible structure interaction”, Comput. Methods Appl. Mech. Engrg.
Vol.198(2009), pp.2650-2661.

(4) Takizawa, K., Christopher, J., Tezduyar T.E. andSathe, S., ”Space-time ﬁnite element computationof arterial ﬂuid-structure interactions with patientspeciﬁc data”, Int. J. Numer. Meth. Biomed. Engng. Vol.26(2010), pp.101-116.

(18) Sugiyama, K., Ii, S., Takeuchi, S., Takagi, S. andMatsumoto Y., ”Full Eulerian simulations of biconcave neo-Hookean particles in a Poiseuille ﬂow”,
Comput. Mech. Vol.46-1(2010), pp.147-157.

(5) Torii, R., Oshima, M., Kobayashi, T., Takagi, K.
and Tezduyar, T.E., ”Fluid-structure interaction
modeling of blood ﬂow and cerebral aneurysm: Signiﬁcance of artery and aneurysm shapes”, Comput. Methods Appl. Mech. Engrg. Vol.198(2009),
pp.3613-3621.

(19) Sugiyama, K., Ii, S., Takeuchi, S., Takagi, S. andMatsumoto Y., ”A full Eulerian ﬁnite diﬀerence
approach for solving ﬂuid-structure coupling problems”, J. Comput. Phys. published online (doi:
10.1016/j.jcp.2010.09.032)
4

c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
B5-4(20) Ii, S., Sugiyama, K., Takeuchi, S., Takagi, S. andMatsumoto, Y., ”An implicit full Eulerian method
for the ﬂuid-structure interaction problem,” International Journal for Numerical Methods in Fluids(DOI: 10.1002/ﬂd.2460)

4

1
0.8

y

(21) Nagano, N., Sugiyama, K., Takeuchi, S., Ii, S., Takagi, S. and Matsumoto Y., ”Full Eulerian ﬁnitediﬀerence simulation of ﬂuid ﬂow in hyperelastic
wavy channel”, Journal of Fluid Science and Technology Vol.5(2010), No.3, pp.475-490.

0.6

2

0.4
0.2

0

Figures

0
0

2

4

6

8

x

Rigid wallg replacements

Fig. 3: Flow ﬁeld at t = 30 for the material of G = 2.
Red: elastic object, white: ﬂuid, blue vectors: velocity
vectors. The velocity vectors are plotted at every 16
and 4 grid points in x and y directions, respectively.

Elastic solidδH0
L(= 2π)ky = H0 (1 − δ sin kx)

Elastic solid

y

h(x)
x

h0

H(x)

H0

10-2

Relative error

Fluid

Rigid wallFig. 1: Schematic of the ﬂow ﬁeld

10-4

Zero-order modeFirst-order mode
Second-order mode
Third-order mode
∝Nx-1
128

256
Nx

512

Fig. 4: Decay of the relatice error with respect to thecase of Nx × Ny = 1024 × 512 for the material of G = 2.

4

y

10-3

2
4

1
0.8

0
2

4

6

8

y

0

x

0.6

2

0.4
0.2

Fig. 2: Initial ﬁber orientation. Grey: elastic object,white: ﬂuid. The ﬁber is plotted on a scale of one to
ﬁve at every 8 and 4 grid points in x and y directions,
respectively.

0

0
0

2

4

6

8

x

Fig. 5: Contour of solid volume fraction ϕs for the caseof G = 10. Red: elastic object, white: ﬂuid. The initial
shape of the elastic wall is indicated by broken line.

5c 2010 by JSFM
Copyright ⃝第 24 回数値流体力学シンポジウム
B5-4

1.4
1.2

y

1
0.8

Simulation results
Theoretical analysis
Initial shape

0.6

0

2

4

6

8

x

Fig. 6: Comparison of the interface geometries between
the simulation results and theoretical analysis for thecase of G = 10 together with the initial geometry of
the wavy elastic wall. Simulation results are plotted atevery 4 grid points in x direction.

1.4

10-2

y

Relative error

1.2

10-3

10-4

Zero-order modeFirst-order mode
Second-order mode
Third-order mode
∝Nx-1
128

1
w/o anisotropy
k1=0.5k1=1.0k1=2.0

0.8
0.6

256
Nx

512

2

4
x

6

8

Fig. 9: Comparison of the interface geometries by varying k1 for the material of G = 2, k2 = 1. Simulation
results are plotted at every 3 grid points in x direction.

Fig. 7: Decay of the relatice error with respect to thecase of Nx × Ny = 1024 × 512 for the material of G = 2,
k1 = 1, k2 = 1.

4

0

1

y

0.8
0.6

2

0.4
0.2

0

0
0

2

4

6

8

x

Fig. 8: Flow ﬁeld at t = 30 for the material of G = 2,
k1 = 1, k2 = 1. Red: elastic object, white: ﬂuid, blue
vectors: velocity vectors, black lines: ﬁber. The velocity
vectors are plotted at every 16 and 4 grid points in x
and y directions, respectively. The ﬁber is plotted on ascale of one to ﬁve at every 8 and 4 grid points in x andy directions, respectively.

6

c 2010 by JSFM
Copyright ⃝