数理解析研究所講究録
第 1808 巻 2012 年 101-120

101

Volume Penalization 法による翼果まわりの流れの数値解析

名古屋大学大学院工学研究科，JST CREST

澤村

陽一

石原

卓

Yoichi Sawamura
Graduate School of Engineering,
Nagoya University
JST, CREST

名古屋大学大学院工学研究科，JST CREST
Takashi Ishihara
Graduate School of Engineering,
Nagoya University

JST, CREST

102
はじめに
植物の中には飛行する種子を持つものがある．例えばタンポポの種子は綿毛を用いて飛行することができる．また，翼を持った種子をもつものもあり，そのような種子を翼果と呼ぶ．カエデの翼果は重心が種子の側に偏っており，木から落下
すると同時に自然に回転する．この回転により，その降下速度は毎秒 1 $m$ と低く抑えられ，カエデは自らの種子を広範囲にひろげることができる．東ら [3](1989)
は，実験および運動量理論による解析の結果，回転しながら落下する翼果は落下速度を抑えるという意味で最適な運動を実現していることを示唆した．また，翼果の質量の大きさに対して翼面積が小さいのにもかかわらず，降下速度が非常に遅いということを明らかにした．

近年，Lentink et al.[1](2009) により，回転しながら落下するカエデやシデの翼
果の翼部分の前縁の上方には，前縁渦 (Leading-Edge Vortices;LEV) と呼ばれる強い渦が現れることが明らかになった．前縁渦により強い揚力が発生し，翼果の降下速度を低下させることに貢献している
また，前縁渦は ${\rm Re}=1000$ 前後のはばたき翼，つまり昆虫のはばたきにおいても発生することが知られている．前縁渦は昆虫の飛行原理のひとつである「失速遅れ」において重要な役割を担っており，前縁渦を解析することは昆虫の飛行原理について理解するために非常に重要である．昆虫のはばたきにおいて発生する前縁渦は発生剥離再捕獲を繰り返す非定常性を持ち，翼果まわりの定常的な前縁渦と比較してその役割はより複雑である．
そこで，本研究では，回転しながら落下する種のまわりの流れにおける翼形状と前縁渦の生成維持機構との関係を数値計算によって明らかにすることを目指
し，フーリエスペクトル法に Volume Penalization 法 [2] を適用した数値計算手法の開発検証と $VP$ 法を用いて行った翼果モデルまわりの流れの数値計算の結果について報告する．

本稿の構成は次のようになっている．第 2 節では，数値計算手法である Volume
Penalization $(VP)$ 法についての説明を行う．第 3 節では，円柱まわりの流れの数$VP$ 法を用い$VP$ 法の検証を行った．第 4 節では，値計算を $VP$ 法を用いて行い，て行った翼果 $2D$ モデルまわりの流れの数値計算の条件および結果について報告する．第 5 節ではまとめを行う．

Volume Penalization 法
翼果のように複雑な物体形状を導入し，かつ翼果の形を変更しての流れの計算$VP$ 法は固体境界を含を容易に行うために Volume Penalization( $VP$ ) 法を用いる．
む流れの数値計算を行う手法のひとつで，埋め込み境界法の一種である．埋め込
み境界法は境界条件を課す代わりに運動方程式の外力項を追加することで物体を表現し，移動変形する複雑形状物体まわりの流れの数値計算に適している．ま103
た，物体適合格子を用意することなく構造格子上で物体を表現することが可能であるため，スペクトル法などの既存の乱流数値計算に用いられた高精度の離散化$VP$ 法は物体を浸透率をもった多孔質媒体とみなす手手法にも適用が可能である．
$VP$ 法は Arquis et al[12]法である．
により，物体上では Darcy の法則が成立し，流体中ではナビエストークス方程式が成立するような手法としてキャビティ内流
れでの数値計算として考案された．その後，Angot ら [5] により一般的な問題に対して拡張された．また，近年は $VP$ 法を含む偽領域法が複数の境界条件をもった条件を含むこと [10], 高次精度を持つこと [11] が証明されている．
法は外力項 (penalization term) を設定し物体形状を表現する手法である．
penalization term を組み込んだ非圧縮 Navier-Storks 方程式は次のようになる．
$VP$

(1)

$\{\begin{array}{l}\partial_{t}u_{\eta}+u_{\eta}\cdot\nabla u_{\eta}+\nabla p_{\eta}-v\nabla^{2}u_{\eta}+X_{\frac{\Omega}{\eta}}(u_{\eta}-u_{s})=f.\nabla\cdot u_{\eta}=0.\end{array}$

ここで， は $VP$ 法を適用した時の速度， は $VP$ 法を適用した時の圧力， は動粘性係数， は物体速度， は外力， は物体の浸透率をそれぞれ示す．本研
は，究では，$f=0$ とする．マスク関数$u_{\eta}$$v$$p_{\eta}$$u_{S}$$f$$\eta$$\chi_{\Omega}$(2)$\chi_{\Omega}=\{\begin{array}{l}1 (X\in\overline{\Omega}_{S})0 (X\in\Omega_{f})\end{array}$

このような式で表現される．マスク関数は，物体中及びその境界上$(\overline{\Omega}_{s})$が 1, 流体となるステップ関数である．
は浸透率を表すパラメータで， がゼロの極限で penalization term を含む方がナビエストークス方程式の解に収束することが Angot らにより解程式の解析的に証明されている [5]. 本研究ではすべての数値計算において $\eta=10^{-3}$ を用い

領域$(\Omega_{f})$が$\eta$$0$$\eta$$u_{\eta}$

る．数値計算手法は流れ関数-渦度法を用いる．流れ関数-渦度法による渦度方程式
は (1) 式の Rotation を取り，次のようになる．

$\partial_{t}\omega_{\eta}+u_{\eta}\cdot\nabla\omega_{\eta}-v\nabla^{2}\omega_{\eta}+\nabla\cross\frac{\chi_{\Omega}}{\eta}(u_{\eta}-u_{s})=\nabla\cross f$.

(3)$\nabla^{2}\psi=-\omega_{\eta}$.

(4)$u_{\eta}=(\partial_{y}\psi, -\partial_{x}\psi)+U_{\infty}$.

(5)

ここで， は $VP$ 法を適用した時の渦度， は流れ関数，$\psi$$\omega_{\eta}$$U_{\infty}$

は平均流速度を示す．空間離散化にはフーリエ・スペクトル法を用い，時間積分には 4 段 4 次精度
Runge-Kutta 法を用いて行う．境界条件は周期境界条件とする．

104
マスク関数の平滑化
フーリエスペクトル法において，マスク関数にステップ関数を用いると Gibbs現象が発生し，数値計算に悪影響を及ぼす．そこで，マスク関数にガウシアンフイ
ルタをかけることでステップ関数を滑らかな形にして計算を行う．ガウシアンフイ
ルタは Kolomenskiy らの研究 [2] を参考に，フーリエ変換を用いて次のようにかける．

$\hat{\chi}=\exp(-C_{smth}((k_{x}^{2}+k_{y}^{2})/N^{2}))\hat{\chi}$(6)

sharp.

方向の格子$N=N_{y}$ とし， はここで， はそれぞれ 方向，y 方向の波数，は2値点数である． は平滑化されたマスク関数をフーリエ変換したもの，$x$$k_{x},k_{y}$$N_{y}$$y$$\hat{\chi}_{sharp}$$\hat{\chi}$

マスク関数をフーリエ変換したものを示す．また，す．本研究では主に $C_{smth}=12.4$ を採用した．

$C_{smth}$は平滑化パラメータを示

物体にかかる流体力の求め方$VP$ 法では応力 $\sigma(u, p)=$は，$(\nabla u+(\nabla u)^{t})-pI$ の計算をすることなく，次の式で求めることができる [5].

物体にかかる流体の力 $F$ およびモーメント$M_{x_{c}}$$F = \int_{\partial\Omega_{8}}\sigma n_{f}d\gamma=\lim_{\etaarrow 0}\int_{\Omega}\frac{\chi_{\Omega}}{\eta}(u_{\eta}-u_{s})d\Omega+V_{c}\ddot{x}_{c}$,(7)$M_{x_{c}} = \int_{\partial\Omega_{S}}(x-x_{c})\cross\sigma n_{f}d\gamma$$= \lim_{\etaarrow 0}\int_{\Omega}\frac{\chi_{\Omega}}{\eta}(x-x_{c})\cross(u_{\eta}-u_{s})d\Omega+I_{c}\ddot{\theta}_{c}$.

(8)ここで， は物体の体積， は断面 2 次モーメント， は物体の重心を示す位
置ベクトル， は物体の回転角をそれぞれ表す．物体の回転と移動については，$I_{c}$$V_{C}$$X_{\mathcal{C}}$$\theta_{c}$Kolomenskiy et al.[2] らによる手法を用いた．この手法はマスク関数のフーリエ変換を行い，その位相シフトを利用することで実行される．

$VP$ 法の検証(円柱まわりの流れの数値計算)法を用いた数値計算が適切に行われていることを確認するため，一様流中に円柱を固定した場合と静止流体中で円柱を移動した場合の円柱の抗力について比$VP$較を行った．

105
数値計算条件

数値計算領域および境界条件を図 1 に示す．また，比較に用いた数値計算条件
を表 1 に示す．円柱直径を , 動粘性係数を , 数値計算領域は 方向，方向そ
$D$

れぞれに

$x$

$\nu$

$y$

方向格子点数を $N_{x},y$ 方向格子点数を, 時間刻みを $dt$ , 一様
流速度を , 円柱移動速度を$C_{smth}$ を 6.2,12.4,24.8,49.6 と変化させ，とする．
$L_{x},$$L_{y},$$x$$N_{y}$$u$$u_{sx}$$t=1.6$ まで数値計算を行った．初期条件では物体・流体ともに静止状態とする．

図 1: 円柱まわりの流れの数値計算の模式図

表 1: 円柱の固定移動時の数値計算条件

数値計算結果
先行研究結果との比較

円柱を固定した場合の抗力と移動させた場合の抗力についての比較を行った結果を図 2 に示す．静止流体中で物体を移動させた場合および，一様流中で物体を$VP$ 法が適切に実装されていると言うこと固定させた場合で抗力が一致しており，
ができる．

106

$\triangleleft oo\alpha$

$0$

0.

$2$

0.

$4$

0.

$6$

0.
time
$8$

1

1.

図 2: 固定円柱と移動円柱の抗力の比較$\triangleleft t9oo$time

図 3: $C$ $mth$ の値による抗力の比較
。

$2$

1.

$4$

1.6

107
抗力の平滑化パラメータ $(C_{smth})$ の変化に対する依存性を 6.2, 12.4, 24.8, 49.6 と変円柱まわりの流れについて，平滑化パラメータ化させた時の抗力の変化を図 3 に示す．平滑化パラメータが減少するにつれ，抗$C_{smth}$

力も減少するという関係がわかった．これは $Re=100$ での無限長円柱の抗力係数から得られる抗力の値に近づく傾向をもつ．
$C_{smth}=12.4$ としたときの抗力は Kolomenskiy らの結果 [2] と定量的に一
また，

致した．

翼果 2 次元モデルの数値計算
Azuma et al の実験 [3] により得られた結果から，翼果のモデルを作成し翼果のまわりの流れの数値計算を行った．東らによる翼果の断面形状，数値計算に用いた翼果 $2D$ モデルの形状を表 2 に示す．表 2 の断面図の列は実物の翼果を上から見た図および赤線の部分でカットした断面のスケッチである．赤線部分は翼果の回転半径の 4 分の 3 の部分 $(r=0.75R)$ を指す．また，測定されたパラメータのうち，レイノルズ数は代表速度を $r=0.75R$ における翼果の回転速度，代表長さを $r=0.75R$ における翼弦長として計算されている．これらのカジカエデ (Acerdiabolicum Blume.) とイロハカエデ (Acer palmatum Thunb.) の断面形状から，
翼果 $2D$ モデルを作成した．翼果 $2D$ モデルは，イロハカエデは 3 点を選び，前縁部を最大翼厚さを直径とする半円で表し，後縁は上三角の直角三角形として表現する．また，カジカエデは前縁部を最大翼厚さを直径とする半円で表し，後縁部を半径 5 の円弧として表現する．

表 2: 観測された翼果の運動と数値計算条件

108(1) 周期境界条件での数値計算
数値計算条件
数値計算領域の模式図を図 4, 数値計算条件を表 3 に示す．数値計算領域は X 方, 一様流速度を
方向格子点数を $N_{x},y$ 方向格子点数を向に方向に$L_{x},$$L_{y},$$y$$N_{y}$$x$$dt=5.0\cross 10^{-4}$ ととする． は翼弦長，が翼厚さ， は迎角を示す．また，する．レイノルズ数は代表速度 を一様流速度 , 代表長さ を翼弦長 として
動粘性係数 を調整することで設定した．

$u=1$

$c$

$\theta$

$t$

$U$

$u$

$L$

$c$

$v$

図 4: 翼果モデルまわりの流れの数値計算領域の模式図

表 3: 数値計算条件

数値計算結果を 5 度から 30 度まで 5 度刻みで変化させたときの揚抗比をプロットしたものである．ここで揚抗比 $LDR$ は次式で定義される．
図 5 は迎角$\theta$$LDR=L/D$

(9)

109
ここで $L$ は揚力で，(7) 式の $F$ の$y$方向成分， は抗力で (7) 式の
$D$

$F$

の $x$ 方向成分をそれぞれ示す．
図 5 より，5 度，10 度の低い迎角では翼上面の渦が剥離せずに残り続け，それに伴い翼の揚抗比も上昇していくことがわかる．一方，15 度よりも高い迎角では最初に揚抗比が上昇するが，時間経過に伴い渦が剥離し，揚抗比が増加したり減少することが分かった．
図 6 はカジカエデのまわりの流れについて，迎角を 5 度から 30 度まで 5 度刻みで変化させたときの揚抗比をプロットしたものである．イロハカエデと異なり，カジカエデでは 10 度の迎角で翼上面の渦が剥離をしはじめることがわかる．また，初期に最大の揚抗比を示す迎角は 10 度である．これは，東らの観測実験で得られた迎角の 12.5 度と近い数値を示している．
イロハカエデ，カジカエデ共に時刻 $t=14$ で迎角によらず揚抗比が急激に低下していることがわかる．これは，一度翼果の後方に流れ去った渦が周期境界条件
によりふたたび翼果の前部から流入し，衝突した時刻に相当している．このことから，周期境界条件下での数値計算に後流が影響を及ぼすことがわかる．

(2) フリンジ領域を用いた数値計算
先述の手法では周期境界条件を課しているため，時間経過により翼後流が前方から流入し流れに悪影響を及ぼすことがわかった．そこで，翼から充分離れた領域で渦を消す外力 (フリンジ) を導入する [4]. 具体的には渦度方程式 (3) の右辺に$\lambda(x)(-\omega_{\eta})$を加える．ここで$\lambda(x)$

(10)

.

は次のように定義される．

$\lambda(x)=\lambda_{\max}[S(\frac{x-x_{start}}{\triangle_{rise}})-S(\frac{x-x_{end}}{\triangle_{fall}}+1)]$

(11)

,

(12)

$S(x)=\{\begin{array}{ll}0 x\leq 01/[1+\exp(\frac{1}{x-1}+\frac{1}{x})] 0<x<11 1\leq x\end{array}$

ここで， はフリンジの強さを示すパラメータである．
fall はフリンジの上昇区間及び加工区間，
はそれぞれフリンジ区間の開始点，終了点を示す． の関数形の例を図 7 に示す．この際のパラメータは $xend-x_{start=}$$\triangle_{r}i_{Se},$$\lambda_{\max}$$x_{start},$$x_{en}d$$\lambda(X)$$4.0,\triangle_{ri\mathcal{S}e=}1.5,\triangle fall=1.0,\lambda_{\max}=10.0$である．

$\triangle$

110$\frac{\underline{o}}{\alpha}$$\approx_{I}\check{-}\alpha^{0}\circ$

5

$0$

10
time

15

20

図 5: イロハカエデのまわりの流れにおける揚抗比 $({\rm Re}=500)$$\check{\alpha}\underline{\circ}$$\check{-4}\triangleleft_{1}\alpha-\omega$

$0$

5

15

$10$

20

time

図 6: カジカエデのまわりの流れにおける揚抗比 $({\rm Re}=1370)$

111

X$\check{d}$$0$0.511.5

22.5

33.5

4

$\cross$

図 7:$\lambda(x)$の例(2-1) フリンジ領域の有無及び翼果に対する位置による影響
フリンジ領域により流れが受ける影響について，迎角 20 度のイロハカエデのまわりの流れでの渦度場及び揚抗比について数値計算を行った．
数値計算条件標を指す．

次の 4 条件で数値計算を行った．ここで，$x_{cg}$は翼果の重心の

$x$

座

条件 A
フリンジ領域を持たない条件 (周期境界条件)
条件 B

フリンジ領域を持ち，フリンジ領域と翼果の位置が近い条件$(フリンジ領域:12 \leq X\leq 16, x_{cg}=8.0)$

条件 C

フリンジ領域を持ち，フリンジ領域と翼果の位置が遠い条件$(フリンジ領域:12 \leq x\leq 16, x_{cg}=3.2)$

条件 D

フリンジ領域を持たず，主流方向に計算領域を拡大し長時間計算を行う条件
条件 A, $B,C,D$ の数値計算条件を表 4 に示す．フリンジ関数についての各パラメー
タは $x_{start=}12,$ $x_{en}d=16,$ $\triangle_{rise=}1.5,$ $\triangle fall=1.0,$ $\lambda_{\max=}10$ のように設定した．

112
条件 A, $B,C$ ではイロハカエデ，カジカエデについて，条件 D ではイロハカエデについてのみ数値計算を行った．

表 4: フリンジ領域の影響を比較した計算条件

数値計算結果 条件 A, $B,C,D$ のそれぞれの渦度場を図 8, 図 9, 図 10, 図 11 にそれぞれ示す．図 8, 図 9, 図 10 は同時刻の渦度場である．図 9, 図 10 を見ると，翼果の後流部分で渦度が消滅し，フリンジ領域が機能していることがわかる．また，翼果の周囲の渦度場を見ると，図 8, 図 10 はほぼ同じ渦構造になっているが，図 9 では後流の渦が潰れたようになっている．このことから，翼近くにあるフリンジ領
域は翼周りの渦構造にも影響を与えることがわかる．条件 D の流れの渦度場では，前縁部と後縁部から渦が交互に剥離して流れ去っていく様子がわかる．
$t=3$
各計算条件における揚抗比の時間変化をプロットしたものを図 12 に示す．のあたりで条件 B の揚抗比が条件 A, と比べて変化していることがわかる．また，$t=4$ 近傍でも同様に揚抗比が条件 A ととで変化していることがわかる．これら
は，フリンジ領域により渦が消滅した影響であると考えられる．このように，フリンジ領域は揚抗比に対して影響を与えることがわかった．また，条件 D が他の条件の揚抗比と異なるのは，領域を拡大したために流れ場の周期性に変化が現れたためであると考えられる．条件 B, $C,D$ ともに揚抗比が振動していることがわかる．この理由は翼の前縁部と後縁部から交互に渦が剥離放出されているためであると考えられる．一方，条件 C, を比較すると揚抗比の時間変化に違いが見られるが揚抗比の変化幅振動数はほぼ同じであることが分かる．以上よりフリンジ
$C$

$C$

$D$

領域を導入し，フリンジ領域から適切に離れた位置に翼を設定することで，フーリエスペクトル法での数値計算を効率的に行うことができるといえる．

113

t $.
$*$

$\mathfrak{W}$

註 59

298
159
189
響6

8

8266

498

688

838

1888

図 8: 条件 A の流れの渦度場$t\# 9.\theta 0$$a\mathfrak{B}{\}$$\mathfrak{W}g$$1\Re$$1\mathfrak{W}$

$\Re$

8296

488

698

888

1988

図 9: 条件 B の流れの渦度場$t\approx \mathfrak{g}.gg$$a\Re$

230
二 58

i33
$\delta\theta$aaos$4\theta\theta$888$\epsilon\epsilon\S$1098

件 の流れの渦度場
図 10: $Ct=$$7Rx$$\lambda 8{\} l\S laew_{0}1l\S \mathfrak{W}$$\mathfrak{g}$$5\theta\theta$$\lambda\theta l\theta$

$A$

$8

$\emptyset$

図 11: 条件 D の流れの渦度場

2 優珀

114

$=\alpha\circ-$

$b0$

邸

$\dot{\underline{\frac{\triangleleft_{1}}{\underline{b}}}}$$0 5 10 15 20$time

図 12: 各条件での揚抗比の比較(2-2) フリンジ領域を含む翼果のまわりの流れの数値計算

これまでの結果を踏まえ，フリンジ領域を含んだ翼果のまわりの流れの数値計算を行った．
(フリンジ領域を持ち，フリ
数値計算条件 フリンジ領域の設定は前節の条件
ンジ領域と翼果の後縁が離れている条件) と同じ条件とする．また，その他の数$C$

値計算条件，翼形状は先に用いた形状と同じとする．数値計算条件を表 5 に示す．

表 5: フリンジ領域を持った翼果のまわりの流れの数値計算条件

数値計算結果 フリンジ領域を含んだ翼果のまわりの流れにおける翼果の揚抗比の時間変動を迎角について比較したものを図 13, 図 14 にそれぞれ示す．

115
図 13 はイロハカエデの揚抗比の時間変動を示す．迎角 5 度，10 度では揚抗比が振動せずに一定になっていることがわかる．また，迎角 15 度では，$t=12$ 以降で揚抗比が振動していることがわかる．迎角 20 度以上では $t=6$ 以降で揚抗比が定常的に振動していることがわかる．渦度場の時間変動をみると，迎角 5 度，10 度では渦が剥離せずに後流へと流れさっていることがわかる．また，迎角 15 度以上では渦が翼の前縁部及び後縁部から周期的に剥離していることがわかる．この剥離が揚抗比の振動に対応していることがわかる．
$t=10$ 以降
図 14 はカジカエデの揚抗比の時間変動を示す．どの迎角においても，で揚抗比が振動しているのがわかる．また，迎角 度を除いて単調な周期の振動を持っていることがわかる．渦度場の時間変動を見ると，これらは翼の前縁部及び後縁部から周期的に渦が剥離することに対応した振動となっている．
図 15 は迎角 10 度のイロハカエデのまわりの流れの渦度場である．図 13 より，
迎角 10 度のイロハカエデのまわりの流れの揚抗比はほぼ一定である．渦度場をみると，翼の後縁から渦が剥離し流れさっている様子がわかる．渦度場は時間変化に対してもほぼ一定であった．
図 16, 図 17 は迎角 15 度のイロハカエデのまわりの流れの渦度場である．図 16は揚抗比が極大の時の時刻，図 17 は揚抗比が極小の時の渦度場である．翼回りの渦構造に大きな違いは見られないが，揚抗比が極小となるときは極大となるとき$1O$と比べて，翼の後縁部から渦が上面にわずかに回り込んでいることがわかる．
カジカエデまわりの流れも含めて観察した結果，次のことが共通した特徴として挙げられる．揚抗比が極大となる場合は，前縁から剥離した渦が後縁付近に再付着をする，あるいは剥離した渦が後縁部上方で大きな渦を形成する場合である．
一方，揚抗比が極小となる場合は翼下面部の渦が翼後縁部より翼上面へと回り込
む場合である．

(2-3) 実験の流れ場との比較
数値計算により得られた準定常状態の渦度場と，Lentink らの実験 [1] により得られた渦度場との比較を行った．比較に利用した数値計算条件は，フリンジ領域を持つ迎角 20 度のイロハカエデのまわりの流れ，及び迎角 10 度のカジカエデのまわりの流れである．比較に利用した数値計算の迎角は，初期の揚抗比が最大となる迎角ではなく，東らにより得られた実験の迎角と近い値を利用した．その結果，イロハカエデ，カジカエデともに揚抗比が極大となる瞬間と実験で定常的に維持されている渦度場には類似性が見られた．一方で，数値計算上では翼面上の渦構造は維持されず，すぐに剥離してしまう．これは流れ場の 2 次元性，つまりス
パン方向の流れ場が存在しないために発生するものであると考えられる．

116$=\circ\infty\infty\circ i)$$-t^{\vee}\underline{\dot{*\triangleleft_{1}.}}$$\deg 05-$$\deg 10 u\cdot\cdot uvx$$\deg 15 \cdots\cdots$$\deg 20/\cdot\cdots\cdots\cdot$$\deg 25$$\deg 30$

$012345678$ 910111213141516171819202122232425
time

図 13: イロハカエデのまわりの流れにおける揚抗比 (フリンジあり)$\check{t9}\underline{o}\zeta\triangleleft bl)$$-\underline{\frac{\mathring {}1}{"}}$$\deg 05-$$\deg 10ru-$$\deg 15$$\cdots\cdots$$\deg 20$$,,,\sim$$\deg 25$$\deg 30$

time

図 14: カジカエデのまわりの流れにおける揚抗比 (フリンジあり)

117$t=22.SB$
$g50$

200

150
100

50
$B$

$\mathfrak{g}$

200

486

606

808

$10B0$

図 15: イロハカエデのまわりの流れの渦度場 $(Re=500, \theta=10, t=22.80)$$t=22.8B$

250
288

158
100
$5B$

$B$

$B$

$2B0$

$4B0$

$BBB$

see

lBBB

図 16: イロハカエデのまわりの流れの渦度場 (揚抗比極大) $(Re=500,$ $\theta=15,$ $t=$
$2280)$$t=24.80$

250

200
158

100
50
$\mathfrak{g}$$p$$g\mathfrak{g}g$$4\emptyset 0$

699

see

$1B09$

図 17: イロハカエデのまわりの流れの渦度場 (揚抗比極小) $(Re=500,$ $\theta=15,$ $t=$

24.80)

118
まとめ
本研究ではフーリエスペクトル法に Volume Penalization( $VP$ ) 法を適用した数値計算手法の検証とその手法を用いた数値計算および解析を実施した．具体的に$VP$ 法の検証を行った．次に，は，まず，円柱まわりの流れの数値計算を実施し，
翼果まわりの 2 次元的な流れの数値計算を実施し，フリンジ領域の導入の影響，および準定常状態における翼果まわりの渦度場と揚抗比の関係を調べ，実験との比較を行った．以下に結果をまとめる．
$VP$ 法の検証$VP$ 法の検証行った．その結果，得ら
円柱まわりの流れの数値計算を行い，れた抗力が先行研究の結果と一致した．また，移動円柱まわりの流れと固定
円柱まわりの流れの抗力が一致した．以上の結果から $VP$ 法が適切に実装されていると確認した．また，平滑化パラメータ $C_{smth}$ を減少させることで円柱まわりの流れの抗力が減少した．これは $C_{smth}$ がゼロの極限で解析解に近づくことを意味している．

翼果 2 次元モデルの数値計算
Azuma et al の実験 [3] により得られた結果から，イロハカエデとカジカエデの翼果のモデルを作成し翼果のまわりの流れの数値計算を行った．
周期境界条件での数値計算

周期境界条件下での数値計算の結果から，次のことがわかった．揚抗比の時間変動は，翼の形状やレイノルズ数の違いにより迎角への依存性が異なることがわかった．イロハカエデでは ${\rm Re}=500$ では
15 度から渦の剥離がはじまり，カジカエデでは ${\rm Re}=1370$ では 10 度から渦の剥離がはじまっていることがわかった．周期境界条件により，後流が数値計
算結果に影響を及ぼし，準定常状態の数値計算が困難であることが分かった．

周期境界条件による後流の影響を排除するために，フリンジ領域を導入した．フリンジ領域を組み込んだ数値計算と，数値計算領域を広げた周期境界条件の数値計算との比較を行い，フリンジ領域を検証しフリンジ領域の導入

た．その結果，フリンジ領域を導入することで翼果のまわりの流れが定常状態になる時刻まで数値計算を行うことができることがわかった．以上から，
周期境界条件を課すスペクトル法を用いた数値計算であっても，フリンジ領域を導入することで翼果の数値計算に利用することが可能であることが示唆された．

フリンジ領域を含む翼果のまわりの流れの数値計算 フリンジ領域を導入した翼果まわりの流れの数値計算により，翼果のまわりの流れは一定のレイノ

119
ルズ数迎角で揚抗比が振動する定常状態となることがわかった．また，揚抗比の振動は渦の剥離に伴い発生することがわかった．揚抗比と渦度場の関係について詳細に調べた．揚抗比が極大となる場合は，前縁から剥離した渦が後縁付近に再付着した場合か，剥離した渦が後縁部上方で大きな渦を形成していた場合であった．揚抗比が極小となる場合は，翼下面部の渦が翼後縁部より翼上面へと回り込み，前縁で剥離した渦の後縁部への付着を妨げた場合であった．揚抗比の極大極小となる場合の渦構造は迎角によらずみられた．
実験との比較

フリンジ領域を導入した数値計算により得られた準定常状態の渦度場と，Lentink らの実験 [1] により得られた渦度場との比較を行った．
その結果，数値計算により得られた準定常状態の揚抗比が最大となる時刻の渦度場と，Lentink らの実験で観測された渦度場が類似した構造を持つことがわかった．一方，実験結果とは異なり，前縁渦は周期的に剥離し，維持されなかった．また，渦の剥離が見られなかった条件にも，先行研究で見られるような前縁部に大きく集中した渦は見られなかった．これは流れ場の 2 次
元性が大きく影響をしていると考えられる．

2 次元の数値計算では渦構造が，実験により得られた渦構造と似た構造を持つことが分かった．一方で，渦構造は維持されず，3 次元の数値計算を行う必要性があることがわかった．また，フリンジ領域を用いた数値計算の結果から，3 次元数値計算においてもフーリエスペクトル法を用いた数値計算が有効であることが示唆された．
今後は 3 次元の数値計算を行い，渦の生成維持機構についての検証を行う．また，遺伝的アルゴリズムなどの最適化アルゴリズムを用いて，翼の形状を変化させながら数値計算を行い，翼果の形状の妥当性についての検討を行いたいと考え
ている．

参考文献
[1] D.Lentink, et al., “Leading-Edge Vortices Elevate Lift of Autorotating Plant
Seeds”,Science 324(2009), 1438.
[2] D.Kolomenskiy, K.Schneider, “A Fourier spectral method for the NavierStokes equations with volume penalization for moving solid obstacles,” Journal of Computatonal Physics, 228 (2009), pp. 5687-5709.
[3] A.Azuma, K.Yasuda, “Flight Performance of Rotary Seeds”,Journal of Theorical Biology 138(1989), pp.23-53.

120[4] LUNDBLADH, A,BERLIN, S.,SKOTE, M., HILDINGS, C. CHOI, J., KIM,
J. & HENNINGSON, D. ., An efficient spectral method for simulation ofincompressible flow over a flat plate.(1999)$S$[5] Angot ., Bruneau C.- ., Fabrie ., “ penalization method to take into
account obstacles in incompressible viscous flows “, Numerische Mathmatik
81(1999) pp. 497-520
$P$

$H$

$P$

$A$

[6] R.B.Srygley& A. L. R. Thomas, “Unconventional lift-generating mechanisms
in free-flying butterflies”,NATURE VOL420(2002),pp. 660-640.

[7] 神部勉，「ながれの事典」，丸善[8] Shizuka Minami,Akira Azuma “Various flying modes of wind-dispersal
seeds”,Journal of Theoretical Biology 225(2003), pp.1-14[9] Yasuda,K., Azuma Akira., “The autorotation boundary in the flight of samaras.”,Journal of Theoretical Biology $185(1997)$ ,pp.313-320[10] A.Sarthou, S. Vincent, J.-P.Caltagirone, P.Angot “Eulerian-Lagrangian gridcoupling and penalty methods for the simulation of multiphase flows interacting with complex objects”, International Journal ofor Numerical Methods
in Fluids 56(8)(2008), pp.1093-1099[11] I. Rami\‘ere, P.Angot, M.Belliard, “ general fictitous domain method withimmersed jumps and multilevel nested structured meshes”, Journal of Computational Physics 225 (2) (2007) pp. 1347-1387$A$[12] E.Arquis,J.-P.Caltagirone, “Sur les conditions hydrodynamiques au voisinage
d’une interface milieu fluide-milieux poreux:application \‘a la convection naturelle”, Comptes Rendus de l‘Academie des Sciences Paris II 299(1984) 1-4.