JHPCN: 学際大規模情報基盤共同利用・共同研究拠点
第13回 シンポジウム ポスター発表 EX22206

直交格子法による移動境界問題の解法に関する研究
佐々木大輔（金沢工業大学），高橋俊（東海大学）

1

目次

1. 移動境界解析：飛翔昆虫の翼解析

2. 固液混相流解析：流体食品の流れ

3. 固気混相流解析：衝撃波による流れ

2目次

1. 移動境界解析：飛翔昆虫の翼解析

2. 固液混相流解析：流体食品の流れ

3. 固気混相流解析：衝撃波による流れ

3

研究背景
𝜔𝑐
𝑘=
2𝑈
k : 無次元周波数 [-]ω : 角速度 [rad/s]
c : 翼弦長 [m]U : 流速 [m/s]https://mintclub.kobe-np.co.jp/skip/event/4089.html➢ コルゲート翼と平板翼の空力特性の比較

4

Corrugated Wing

一様流・翼弦長基準Re数2800

無次元周波数0.1～0.9

羽ばたき運動ヒービング（上下方向のみ）

長沼，佐々木，高橋，岡本，
第53 回流体力学講演会／第 39 回航空数値シミュレーション技術シンポジウム
講演論文集，2021年

計算手法
支配方程式
計算格子
物体表現
物体境界条件5

非圧縮性3次元NS方程式
等間隔直交格子レベルセット法
ゴーストセル法

計算結果（上昇時の渦度）
コルゲート翼
k=0.2

6

平板

コルゲート翼k=0.4

平板

計算結果（上昇時のCp）
コルゲート翼
k=0.2

7

平板

コルゲート翼k=0.4

平板

計算結果（揚力係数）
0.6
0.4

0.6

k=0.2

0.2

0.2

0Cl

Cl

k=0.4

0.4

0-0.2

-0.2

-0.4

-0.4

-0.6

-10-505

10effective AOA [deg]flat plate (1.5%c)-0.6-10-50effective AOA[deg]5corrugated wing

8 ➢ いずれの場合にもコルゲート翼の方が高Cl

10

目次

1. 移動境界解析：飛翔昆虫の翼解析

2. 固液混相流解析：流体食品の流れ

3. 固気混相流解析：衝撃波による流れ

9

解析手法
支配方程式

周期境界条件

連続の式

D𝛻∙𝑈 =0r
3次元非圧縮性Navier-Stokes方程式
𝜕𝑈1
𝜇 2+ 𝑈・𝛻𝑈 = − 𝛻𝑃 + 𝛻 𝑈𝜕𝑡𝜌𝜌

Ptop= Pbottom -ΔP
流速U
壁面:すべり無しzL
格子

完全等間隔直交格子

離散化

有限差分法

流体-物体境界条件

埋め込み境界法
Level set法，Ghost cell法

圧力解法

BiCGSTAB法

粘性項

2次精度中心差分法

対流項

5次精度WENO法

10

粘性項の時間積分

Crank-Nicolson法

対流項の時間積分

3次精度TVD Runge-Kutta法

yw
UxPbottom
流速U初期条件:Hagen–Poiseuille流れ2
速度分布定義𝛥𝑃 𝐷𝑤(𝑟) =− 𝑟24𝜇𝐿 2
周期境界条件：底面圧力、上面圧力の差一定
流入口，流出口速度が周期境界条件

埋め込み境界法
複雑形状･移動物体定義に有効
流体領域： Fluid Cell（FC）
固体： Object Cell（OC)
流体と固体の境界： Ghost Cell（GC)

埋め込み境界法:
格子上のlevel set法 ,ghost cell法
物体,流体境界定義
cellの分類
𝑑𝐹𝐶 > 0𝑑𝐺𝐶 ≤ 0 𝑎𝑛𝑑 𝑑𝐺𝐶 ≥ −2.25∆𝑥𝑑𝑂𝐶 < −2.25∆𝑥

ghost cellに付加される速度11𝑑𝐼𝑃 + 𝑑𝐺𝐶
𝑈𝐺𝐶 = 𝑈𝐼𝑃 −(𝑈𝐼𝑃 − 𝑈𝐼𝐵 )𝑑𝐼𝑃cellの分類

埋め込み境界法

S. Takahashi, T. Nonomura, and K. Fukuda, Journalof Applied Mathematics, 2014 (2014),pp.1-21.

円管の可視化(α=0.15)
球分割数

格子幅[m]

格子点数

200.15×10-3

100×100×466=4,660,000

12混相流の定量比較(α=0.15)
固体占有率:α球分割数

格子幅[m]

格子点数

100.30×10-357×57×233=757,017

120.25×10-365×65×280=1,183,000

15

0.20×10-375×75×350=1,968,750

200.15×10-3

100×100×466=4,660,000

0.28m[kg/s]

0.26

0.24

0.22

130.15[mm]

実験値とCFDの良好な一致を確認

0.05

0.1

0.15 0.2 0.25
gridsize[mm]0.3質量流量と格子幅の関係(α=0.15)

0.35

目次

1. 移動境界解析：飛翔昆虫の翼解析

2. 固液混相流解析：流体食品の流れ

3. 固気混相流解析：衝撃波による流れ

14

支配方程式

continuous eq.

𝜕𝜌 𝜕 𝜌𝑢𝑖+=0𝜕𝑡𝜕𝑥𝑖

Sutherland’s law

momentum eq.

𝜕 𝜌𝑢𝑖 𝑢𝑗𝜕 𝜌𝑢𝑖
𝜕𝑝 𝑀𝑎 𝜕 𝜏𝑖𝑗+=−+𝜕𝑡𝜕𝑥𝑗
𝜕𝑥𝑖 𝑅𝑒 𝜕𝑥𝑗

𝑇෨
𝜇=𝑇𝑟𝑒𝑓32𝑇
𝑟𝑒𝑓 + 𝑆𝑇෨ + 𝑆

energy eq.

𝜕 𝜌𝐸𝜕𝑀𝑎 𝜕 𝑢𝑖 𝜏𝑖𝑗𝜕𝑞𝑖+𝜌𝐸 + 𝑝 𝑢𝑗 =−𝜕𝑡𝜕𝑥𝑗
𝑅𝑒 𝜕𝑥𝑗𝜕𝑥𝑖
equation of state𝑝= 𝛾−11
𝜌𝐸 − 𝜌𝑢𝑖 𝑢𝑖
2stress tensor

𝜕𝑢𝑖 𝜕𝑢𝑗 2 𝜕𝑢𝑘
𝜏𝑖𝑗 = 𝜇+−𝛿𝜕𝑥𝑗 𝜕𝑥𝑖 3 𝜕𝑥𝑘 𝑖𝑗

heat flux

15𝛾 𝑀𝑎 𝜕 𝑝𝑞𝑖 = −𝜇𝛾 − 1 𝑅𝑒𝑃𝑟 𝜕𝑥𝑖 𝜌𝑝𝑇=𝛾
𝜌equation of motion𝜌𝜋𝑠 𝜕𝑈𝑗=න
6 𝜕𝑡Γ

𝜕𝑢𝑖 𝜕𝑢𝑗−𝑝𝛿𝑖𝑗 + 𝜇+𝜕𝑥𝑗 𝜕𝑥𝑖𝑛𝑗 𝑑𝑆• 無次元化Navier-Stokes方程式• Sutherland則• 球の3自由度並進の運動方程式

計算手法
流体の方程式

圧縮性Navier-Stokes方程式 & Sutherland式

運動の方程式

3自由度並進運動方程式（回転含まない）

空間離散化

等間隔直交格子（外部境界近傍は不等間隔直交格子）

非粘性流束

3次精度MUSCL-Roe & 2次精度擬混合型ハイブリッド

粘性流束

2次精度中心差分

時間積分

3次精度TVD Runge-Kutta陽解法

壁面境界条件ゴーストセル型埋め込み境界法（断熱・粘性壁）

外部境界条件

流入：Dirichlet型 側面&流出：Neumann型

16

ハイブリッドMPI

並列計算法

埋め込み境界法（IBM）
Takahashi, S., et al.,J. of Appl. Math. (2014) ,252478

Nonomura, T., et al.,Math. Prob. Eng. (2017) ,8314615

Mizuno, Y., et al.,Appl. Sci. (2018) ,8, 2387

➢ ゴーストセル型IBMの開発➢ 表面流体力計算の開発➢ 運動方程式連成の開発

17亜音速～超音速の低Re数の物体周りの数値計算法の開発

先行研究に基づいた初期条件
Tanno, H., et al.,Shock Waves (2003) 13: 191–200➢ 衝撃波Ma数：1.22➢ 粒子Re数：30万➢ 衝撃波管による実験&CFD➢ 論文の条件抜粋

Sun, M., et al.,Shock Waves (2004) 14: 3–9➢ 衝撃波Ma数：1.22➢ 粒子Re数：49~49万➢ BFC軸対称細分化格子CFD➢ 論文の記述を元に無次元化

18高圧領域(2)：窒素低圧領域(1)：室温空気

ρ2 = 1.66 kg/m3
u2 = 114 m/s
p2 = 159 kPa
μ2 = ?

ρ1 = 1.205 kg/m3
u1 = 0D=80mm p = 101 kPa
1Exp
μ1 = 1.82x10-5 Pa·s

高圧領域(2)：理想気体

低圧領域(1)：理想気体

ρ2 = 1.9269
u2 = 0.3361
D=1p2 = 1.5696
μ2 = 0.013217 Sun≒本研究ρ1 = 1.4
u1 = 0p1 = 1
μ1 = μ2

本研究は粒子Re数が49の場合のみ(Kn=9.4×10-3)計算格子※CPD: Cell Per Diameter※1/8粗視化データを可視化

単一球の解析

名称

格子幅

CPD20 0.05D

19

格子点数
1.3M (160×90×90)

2球の解析

名称

格子幅

格子点数

球の比重0.025D(CPD40)26M(460×240×240)

10, 100, 1000
(気体密度に対する比)直列

CPD40 0.025D

5.1M (260×140×140)

並列

CPD80 0.0125D

26M (460×240×240)

斜め【検証】格子依存性

圧力の可視化

t=2.0
CPD020

CPD040

CPD080

CPD020Sun

CPD040

CPD080

Sun

全抗力
CPD:
Cell Per Diameter

646420020
246nondimensional time810

CPD080

Sun

摩擦抗力

10820

圧力抗力
friction drag coefficient

8

CPD0401210

pressure drag coefficient

10

t=6.0

CPD0201212drag coefficient

t=4.0864200246nondimensional time

8100246nondimensional time

CPD40で格子収束→移動球の解析はすべてCPD40で実施

810※ここでの比重→気体に対する密度比【応用2】比重別の移動する２球の解析

比重10

比重100021

比重1000は比重10に比べて固定と見做して差し支えない（移動量100倍差）

衝撃波負荷により移動する多数粒子

22

3粒子

16粒子

まとめと今後の課題

1. ゴーストセル法の応用性を移動境界問題で確認

2. ゴーストフルード法を非圧縮・圧縮流れに適用
3. さらに高精度かつ汎用性の高い解析法の構築

23