Kobe University Repository

Kobe University Repository : Thesis
学位論文題目
Title
埋込境界法の一般化による流れの数値シミュレーション
法への応用
氏名
Author
千秋, 雅信
専攻分野
Degree
博士（工学）
学位授与の日付
Date of Degree
2012-03-25
資源タイプ
Resource Type
Thesis or Dissertation / 学位論文
報告番号
Report Number
甲5493
URL
http://www.lib.kobe-u.ac.jp/handle_kernel/D1005493
※当コンテンツは神戸大学の学術成果です。無断複製・不正使用等を禁じます。
著作権法で認められている範囲内で、適切にご利用ください。
Create Date: 2017-03-30
博士論文
埋込境界法の一般化による流れの
数値シミュレーション法への応用
平成 2
4年 1月
神戸大学大学院工学研究科
千秋雅信
論文要旨
日本は山がちな独特の地形を有しており，構造物に対する耐風設計はもちろ
ん風力発電の湿地候補地選定にも複雑な地形と人工構造物の影響を考慮した風
況の予測が必要になる.
複雑な形状の境界を持つ流れ場を数値解析する場合，地形に沿うように形成
される境界適合格子を用いる方法や，複数の異なるメッシュを用いる重合格子
法，物体境界の貫通により生じる格子の破片上で計算を行うカットセル法や物
理変数を境界上で定義し，計算点から境界までの距離を離散化の際に用いる
S
h
o
r
t
l
e
y
-W
e
l
l
e
r近似等の手法があり，それぞれコーディングの煩雑さ，計算の効
率などの点で一長一短である.
近年，複雑な格子形成を行わず運動方程式に外力を導入し境界を表現する
IBM (
ImmersedBoundaryMethod) が種々の流れの数値解析に応用されその
有用性が示されている.元々は心臓内の血流解析や血管内流れなど移動境界を
扱う際に用いられていたが，その後固定境界への応用がなされ最近では L
ESを
用いて丘地形を超える流体解析に用いられ，その際境界条件に対数則を導入し
たものも見られる.この方法は直交格子でしかも粘着条件など境界で値をゼロ
に設定する場合演算が非常に簡単になる.従ってこれまでの IBMの応用の殆ど
tal.は格子で解像困難な
は，直交座標で速度ゼロの境界に限られている. Cuie
形状に IBMを適用している.
丘陵地形とその上の構造物はともに境界形状を複雑にし，その回りの気流性
状に影響を与える.大きなスケールで穏やかに変化する地形は緩やかではある
が増速域を上空にまで広げるなど広域に変化を及ぼす.それに対し局所的に突
出た構造物は，限られた領域ではあるが，逆流や変動など大きな変化を起こし，
それが周囲に伝わる.前者のような地形上の流れは従来から地表摩擦速度や粗
度高さをパラメータとした相似則あるいはそれからのずれの修正で巧く表現さ
れる.従ってその影響は境界(地表)に沿う座標を用いることで的確にまた有
効に考慮できる.しかし，それらには地表面からの距離の対数関数や，陰的に
決定される粗度高さや粗度関数など未知数である速度の非線形関数が用いられ，
速度を指定する IBMでは対応は非常に複雑になる.後者の小スケール物体であ
るが，曲線座標をこれら角ばった形状に沿わせることは困難であるのみならず
不可能な場合があり，またこれらが格子で適合できたとしても物体表面上の境
界条件はより大きなスケールで成り立つ壁面則ではなく，物体位置で流れが遮
蔽される速度ゼ口条件である.従ってこういった構造物等物体の影響は IBMで
の対応が適している.
1
また、この IBMの利点を運動方程式以外の輸送方程式にも応用できれば、熱
や物質拡散の計算にも使えると期待される。本研究では温度や濃度など速度場
だけでなく Nuemann条件やその他の条件に適用できる IBMを提案し、計算例
を示すことによりその適用性、実用性などを検討した。この方法では一般化し
た IBMは乱流エネルギーや散逸などの乱流統計量の計算にも応用可能である。
本研究では従来から風況予測の際に用いられてきた境界適合曲線座標と直交
格子を用いて複雑な境界流れに応用されてきた IBMを融合させ，地形上に存在
する構造物を考慮した気流解析手法の開発と計算手法の検証・応用を行い，以
下のような結論を得た.
①境界適合曲線座標に IBMを融合させたプログラムの妥当性を示すために地
形上の任意の位置に存在する任意の大きさの 2次元・ 3次元形状の構造物
境界周りの流れを矩形近似法と比較し，その妥当性を得た.
②実地形上に存在する建物周りの流れの解析を行い，従来法では不可能だ、つ
た境界表現が可能になり，境界適合曲線座標を用いても IBMの境界表現が
有用であるという見識を得た.
実現象をとらえる場合，その流れはレイノルズ数の高い乱流である.その際，
境界面においては粗度を考慮、したり，対数則を用いなくては実現象をとらえる
ことは不可能である.しかし構造物・建物などの境界は流れが直角にぶつかり，
流れが止まることを考えるとノンスリップ条件を用いることが妥当だと考えら
れる.従来の方法では矩形近似法により境界を表現する手法が見られた.矩形
近似法は直交格子を用いることで計算プログラムはとても簡単なもので済むが，
地表面においてもノンスリップ境界しか与えることができない.また，高い Re
数を対象とする大気流の地表境界条件に IBMを用いることも同様の理由で適し
ているとは言えない.
本研究の計算手法は，地上に存在する建物の境界はノンスリップ条件を IBM
で表現し，地表面の境界については境界適合曲線座標で表現しているため、地
表面の境界条件に粗度を与えたりすることにより対数則を導入することも可能
であり、今後の風況予測において実用的なものとなり得る.
また、通常速度についての運動方程式に用いられる IBMを一般化し，温度場，
混合物，濃度場の境界条件にも適用できるものを提案した.この一般化した IBM
を浮力壁噴流の LES計算に適用し，実験結果及び k
ε モデルを用いた RANS
計算結果と比較しその妥当性を確認した上で，断熱物体，等温物体のある場合
の流れ場と温度濃度場の LES計算を行った. LESのサブグリッドモデル，境界
の取り扱い，さらには解像度等の詳細を検討する必要はあるが，計算結果は概
ね良好であり，計算負荷の小さいことも判明した.今後は，河川や沿岸域での
温水や下水処理水の放流や拡散の予測への活用も期待できる.
1
1
目次
l.序論.•
2
. 埋込境界法について...................................................... 4
2
.
1f
e
e
d
b
a
c
kf
o
r
c
i
n
g法.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 5
2
.
2D
i
r
e
c
tf
o
r
c
i
n
g法.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 6
2
.2体積力の内挿方法.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 7
(1)線形内挿.• • • • • • • • • • • • • • • • .• • • • • • • .• • • • • • • • • • • • • • • • • • • .• • • • • • • • • • • • •• 8
(2)双一次内挿.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 9
(
3
)三次元場での内挿.• • • • • • • • • • • .• • • • • • • • • • • • • • • • • • • • .• • • • • • .• • • • • •
1
0
3
. 境界埋込法の 般曲線座標系への適応.• • • • • • • • • • • • • • • • • • • • .• • • • • • .• • • • • • •• 1
2
3
.1一般曲線座標での基礎方程式 .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
1
2
3
.2シグマ座標系における微分係数の取り扱い.• • • • • • • • • • .• • • • • .• • • • • • • • • • • •• 1
2
3
.3空間差分法.• • • • • • • • • • • • • • • .• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 1
3
4
3
.
4時間進行法.• • • • • • • • • • • • • • • • ..• • • • • • • • • • • • • • ..• • • • • • .• • • • • • • • • • • • • .• •• 1
(
1
)C
r
a
n
k
N
i
c
o
l
s
o
n法.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 1
4
(
2
)A
d
a
m
s
B
a
s
h
f
o
r
t
h法.................................................
. 1
5
3
.5F
r
a
c
t
i
o
n
a
lS
t
e
pM
e
t
h
o
dにおける埋込境界法の導入.• • • • • • • • • • • • • • • • • • • • • • •• 1
5
8
3
.
6計算格子.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • .• • • • • • • • • • • • • • • • • • • • • • • • • • •• 1
4
9
4
. 数値計算手法の検証.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • .• • • • • • .• • • •• 1
4
.1計算条件.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 1
9
4
.
2計算結果と考察.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 2
1
B
Mによる 2次元障害物のある地形上気流解析.• • • • • • • • • • • • • • • • •• 2
6
5
. 一般座標系 I
5
.1計算対象.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • .• • • • • • • .• • • • • •• 2
6
5
.2計算結果と考察.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • .• • • • • • • • • • • • • • • • • • • •• 2
7
(1)境界と速度定義点が一致している場合 (
c
a
s
e1
) ••••••••••••••••••.••••• 2
8
(
2
)境界と速度定義点が一致していない場合 1(
c
a
s
e
2
).
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.3
0
(
3
)境界と速度定義点が一致していない場合 2(
c
a
s
e
3
).
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.3
3
6
. 一般座標系 I
B
Mによる 3次元構造物のある地形上気流解析.• • • • • • • • • • • • • • • • •• 3
6
6
.1計算対象.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • .• • • • • • • .• • • • • • • • • • • • • • •• 3
6
I
6
.2計算結果と考察.• • • • • • • • • • • • • • • .• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 3
6
(1)境界と速度定義点が一致している場合 (
c
a
s
e1
) ••..•••••••••••••••••••• 3
6
c
a
s
e
2
).
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.4
0
(
2
)境界と速度定義点がずれている場合 (
c
a
s
e
3
).
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.4
3
(
3
)構造物が複数存在する場合 (
7
. 一般座標系 I
B
Mの実地形上気流解析への適用.• • • • • • • • • • • • • • • •• • • • •• • • • • • • •• 4
7
7
.1計算対象.• • • • • • • • • • • • • • • • • • • • • • • • .• • • • • • • • .• • • • • • • • • • • • •• • • • •• • •• • • •• 4
7
7
.2計算結果と考察.• • • .• • • • • • • • • • • • • .• • .
.
.• • • • •• • • .• • • • • • • • • •• • • • • • • •• ••• 4
8
8
. 境界埋込法の一般化.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 5
5
5
8
.1N
e
u
m
a
n
n型境界条件への適用..• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 5
8
.2運動方程式以外の輸送方程式への応用.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 5
6
8
.3浮力噴流の L
E
Sシミュレーション.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 5
7
(1)平坦底面上の壁噴流の計算.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• 5
8
(
2
) 平坦底面に断熱障害物のある場.• • • • • • • • • .• .• • • • • • • • • • • • • .• • • • • • • • • • •• 6
9
(
3
)平坦底面に定温断障害物がある場合.• • • • • .• .• .• .• • • • • • • • • • • • • • • • • • • • • •• 7
6
9
. 結論.• • • • • • • • • • • • • • • .• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • .• • • • • •
8
3
謝辞........................................• .........................• ..• 8
5
参考文献.• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • .• • • • .• • .• • • • • • • • • • • • • • • • • • • • •• 8
5
I
I
第 1章 序 論
高度経済成長期に大量じ建設された社会基盤となる構造物が溢れかえる世の中で，近
年台風や突風といった風による災害や、ゲリラ豪雨等の局所的な集中豪雨が相次ぎ，去
1 日に発生した東日本大震災による巨大津波の被害等様々な災害が日
る平成 23年 3月 1
本の国土や人々の生活を脅かしている.一方、東日本大震災の福島原発事故の影響もあ
り、これまで以上に「持続可能な社会の構築」が注目されており，環境に配慮したエネ
ルギーとしての風力発電，都市部における環境汚染物質の拡散など，人々が安心安全に
暮らしを考えていく上で大気環境の適切な予測は不可欠であり，多方面において必要と
されている.
大気環境の予測といえども様々なものがある.大スケールでは，水平距離数百 krnに
r
nから数百 m の領域を扱う局所スケールの流れまで様々であ
も及ぶ気象モデルから数 k
る.気象モデルのような超広域を扱う際，構造物などの流れに対する障害物の影響は
微々たるものであるが，局所スケールの領域を扱う場合は，地形上に存在する樹木や構
造物の流れに対する影響も考慮せざるを得なくなる.風況予測法としては， MASCOT1)
の開発により全国任意地点の詳細な風況の予測が可能となり，その予測誤差も地域気象
モデルにより従来の方法の 3分の l程度まで減少するなどの成果を収めている.この
MASCOT は非線形風況予測モデルであり，一般座標系を用いることで任意傾斜角度の
地形に応用できるものである.海外においても WASp2)や ARPS3)などの局所風況予測法
により都市部の大気の流れ解析 4)に利用されたりしている.しかし，海外と比べて地形
が複雑な日本においては，平坦な地形あるいは緩やかな斜面程度でしか精度が確保でき
ないなどの難点がある. NEDOの研究プロジェクトの一環として開発された LAWEPS5)
は急峻な山岳地帯の風況予測においても年平均風速を 10%以内の誤差にとどめること
に成功している.
風況予測に用いられる方法としては大きく三つあり，一つ目は観測である.風速や風
向を現地で観測することにより現実の流れ場を把握することができる.しかし，観測に
より一度に得られるデータは限られており，現実の流れを完壁に把握することは多大な
コストに加え，長期間の断続的な観測が必要となる.二つ目は風洞実験である.この手
法は長年にわたり利用され続けており，現在ではその精度や信頼性は向上している.し
かし，詳細な模型が必要となる実験は時間と費用がかかり，風洞の大きさや，使用機材
によってもその精度や再現規模が制限されてしまう.三つ目は数値計算である.計算機
器の飛躍的な発達もあり，より広範囲の流れ場を扱うことも可能となってきている.ま
た，実験と比べて容易に条件(領域に流入する流体の速度やその向き等)を変更するこ
とが可能であり，一度の計算で得られるデータ量も多い.このように，それぞれの手法
に一長一短があり，計算の検証に風洞実験や観測の結果が用いられもする.中でも，特
に数値計算はコンビュータの高速化，計算方法の改善などによりその可能性は高まって
きており，風況予測の際の数値計算の重要性も高まっている.
i
m
u
l
a
t
i
o
n
)法
， RANS(Reynolds-averaged
解 析 手 法 と し て は DNS(Direct Numerical S
， LES(Large-EddyS
i
m
u
l
a
t
i
o
n
)法があり， DNSは基礎式に
N
a
v
i
e
r
S
t
o
k
e
se
q
u
a
t
i
o
nm
o
d
e
l
)法
よって与えられる解を，付加的な数値モデルを用いることなく忠実に再現するものであ
る.だが，この方法は解析の際に流れにおける最小スケールまで分解しなければならず，
膨大な格子数と細かな時間間隔が必要となり本解析のような複雑地形上の気流解析に
は適さない. RANS 法は工学的な諸問題に対応し得る最も実践的な手法とされている.
この手法は基礎式を時間平均化したレイノルズ方程式のレイノルズ応力項をモデルに
より与えて数値解を求めるが，解析のほとんど、が経験モデ、ルに頼っており普遍性も低い.
そして本解析で用いた LES 法は，基礎式を格子平均化し，格子幅より大きいスケール
の変動は直接，格子幅より小さいスケールの変動はモデル化することで与える手法であ
る.工学的な流れ場における高レイノルズ数流れに対しては有効とされている.
日本は国土が山がちであり，世界でも独特の地形を有しており，構造物に対する耐風
設計はもちろん，風力発電の風車の設置候補地選定の際にも，風況解析の際に複雑な地
形をどのように扱うかという問題を抱えることとなる.風況予測の際には自然地形の山
や丘，谷などによる影響だけではなく地形上に存在する構造物による風況への影響も捉
える必要がある場合も多々存在してくるであろう.
複雑な境界を扱う方法としては，地形に沿うように形成された格子(境界適合格子)
を用いる方法があげられるの.しかし，境界適合格子では境界地形が複雑になればなる
ほど格子生成が困難になり，地形上に存在する人工的な形状を有する構造物を考慮する
)を用いれば
様な格子形成は不可能である.複数の異なるメッシュを用いる重合格子法 7
地形上に存在する人工的な形状を有する構造物を考慮することは可能だが，異なる格子
の結合部での物理量の補完がかえって煩雑になってしまうなどの問題点があり構造物
が多数存在する領域を対象にするのは現実的で、はない.物体境界の貫通により生じる格
子の破片上で計算を行うカットセル法 8)や物理変数を境界上で定義し，計算点から境界
h
o
r
t
l
e
y
-W
e
l
l
e
r近似 9)においても，前者には極端に小
までの距離を離散化の際に用いる S
さなカットセルの出現や，カットセルと隣接格子との内挿などの課題が.後者にも離散
化が複雑になるなど、の課題が残ってしまう.さらに，カットセル法， S
h
o
r
t
l
e
y
-W
e
l
l
e
r近
似は主に直交格子を用いて複雑な境界を扱う方法として開発されており，境界適合格子
との融合はきわめて困難であると考えられる.
複雑な格子形成を行わずに N
a
v
i
e
r
S
t
o
k
e
s方程式の外力項により境界を表現する埋込
境界法 CIBM(
ImmersedBoundaryMethod)) も種々の流れに適応されている
心臓内の血流解析
1
0
)
元々は
1
1
)や血管内流れ 1
2
)など移動境界を扱う際に用いられていたが，ぞの
後固定境界への応用がなされ，最近では LES を用いて丘地形を超える流体解析にも用
いられ，その際境界条件に対数則を導入したものも見られる
1
3
)
これらの背景を踏まえて本研究では，境界適合格子と埋込境界法を融合して用いるこ
とで，効率よく正確に地形上に存在する構造物・建物の境界を表現する手法を構築した.
この手法を用いて複雑地形上において解析を行うことで手法の妥当性を示し，実地形上
での気流解析に応用した.
この埋込境界法の利点を運動方程式以外の熱や物質拡散の輸送方程式にも応用でき
れば，河川や沿岸への温水や混合物質の拡散の予測法の改良につながると期待される.
そこで本研究では，上述の速度場だけでなく，温度や濃度など Neumann条件やその他
の条件に適用出来る埋込境界法を提案し，計算例を示すことによりその適用性，実用性
などを検討する.この方法で一般化した埋込境界法は乱流エネルギーや散逸などの乱流
統計量の計算にも応用可能である.
本論文は 8つの章により構成されている.第 l章では，研究背景と本研究の目的を述
べた.第 2章では，境界表現手法である埋込境界法ついて，そして任意位置の境界を扱
う際の内挿方法について述べる.第 3章では，計算に用いる基礎方程式，計算アルゴリ
ズムなどについて述べる. 第 4章では，本研究に用いる計算プログラムの妥当性を， 2
次元丘地形を越える低レイノルズ数流れの計算により検証する.第 5章では
2次元
丘地形の上に 2次元形状の構造物を設置し，それらを越える計算を矩形近似法と比較し，
境界での平均流速分布などの考察を行うことで，本研究に用いた境界表現手法の妥当性
について検証する.第 6章では
2次元丘地形の上に 3次元形状の構造物を設置し，
それらを越える計算を矩形近似法と比較し，境界での平均流速分布などの考察を行うこ
とで，本研究に用いた境界表現手法の妥当性について検証する.第 7章では，実地形上
での気流解析に応用し，課題などを検討する.第 8章では，埋込境界法を拡張し，速度
がゼ口になる粘着条件
e
r
i
c
h
l
e
t 条件ではなく，法線方向勾配を指定する
いわゆる D
Neumann 条件への適用法を示す.これは熱や物質拡散の地表境界条件に相当し本方法
が速度についての境界条件のみでなく他の輸送式にも使えることを意味している.さら
に，壁面近傍で噴射される熱噴流に適用し，温度一定の壁や断熱壁の条件が効率よく設
定出来ることを示す.第 9章は，本研究のまとめとしての考察を行い，結論とする.
第 2章 埋 込 境 界 法 に つ い て
複雑な形状の境界を有する流れの解析を行う際，その境界の扱いが課題となる.境界
適合格子などを用いても，ある程度の複雑な形状は扱うことができるが，自然地形に近
い形の境界を扱うまで煩雑化するとこの手法の限界が見えてくる.重合格子を用いるに
しても，いくつもの異なる形状の格子を扱う際にはかえって異なる格子の結合部での物
理量の補完が煩雑となりどちらも限界がある.埋込境界法は，そういった複雑な境界の
扱いを直交格子を用いて可能にした手法である.
埋込境界法は決して新しいものではなく、 40年程前から多くの研究者がこの手法を用
e
s
k
i
nが低いレイノルズ数に仮定された二次元心臓内の僧帽
いている.1970年代初めに P
弁における血流の再現に試みている.その後 McQueenと PeskinI4)15)により三次元手法に
拡張され，心臓や収縮境界の再現に用いられた.上述の研究では，境界の挙動は流れと
既知である境界の聞に働く力のみで決定されることが特徴である.その後も G
o
l
d
s
t
e
i
n
などにより円筒形の構造物周辺流れの解析に用いられるなどしたが、流れと境界の聞に
働かせる力(体積力)の表現に自由定数が存在し，その与え方によっては方程式を解く
際に手間がかかる事もあるようである.しかし，境界適合格子法などに比べると非常に
使い勝手が良いため，現在でも三次元計算に使われている.
その後，体積力の表現を改良することにより計算の安定性に影響しない手法が
Mohd-Yusofにより開発された
1
6
)
G
o
l
d
s
t
e
i
n らの方法で存在していた自由定数が存在し
ない手法であり，効率性の面ではこちらの方法が優れており，三次元複雑流れへの適応
性もあるといえる.これまでにも，曲線形のノズルから出る渦環構造や球周辺の三次元
線対称流れ，モーターを備えた ICピストン内の三次元乱流などに応用されている.
t
)を時刻 tで境界位置
埋込境界法の原理は，位置と時間の関数である体積力点Xb，
ま
Xb
たは境界内部に作用させる事により境界位置で境界条件を満足させることである.体積
力fは流速によるので非定常であり，境界が移動境界の場合は Xbもfも時間の関数とな
る.体積力 fを導入した N-S運動方程式は
竺+u.Vu=-Vp+vV2U +f
d
t
(
2・1
)
である.
この式は fが既知であれば通常の N-S運動方程式に準じた方法で解く事ができる.図
2
.1は計算格子と境界の位置関係を表している.セル中央が速度定義点であり，赤い線
2
.1 (
a
)は境界と速度の定義点が一致
が境界を，紫の領域が固体領域を表している.図 4
している場合であり，直接その点において体積力を作用させて境界条件 u=O(
n
o
n
s
l
i
pの
時)を満たすことができる.しかし，図 2
.1(
b
)のように境界と速度の定義点が一致し
ていない場合は境界内部第一点目に境界位置で境界条件を満たすような体積力を作用
2
.1(
b
)の場合では，流体領域の U
1と境界内第一点目の U
2を内挿
させることとなる.図 して境界で u=Oを満たすように U
2を発生させる体積力を決定しなくてはならないという
ことになる.
UJ
-惨
U
、
r
圃
‘
世 間 伽c
i
r
ゆ
2
(伽n
f
i
町均)
一
園
一
(
a
)境界と速度定義点が一致
(
b
)境界と速度定義点が不一致
2
.1 埋込境界法法概略
図e
e
d
b
a
c
kf
o
r
c
i
n
g法と D
i
r
e
c
tf
o
r
c
i
n
g
先に述べたとおり，体積力の与え方は二通りあり， f
法とがある.
2
.
1 f
e
e
d
b
a
c
kf
o
r
c
i
n
g法
G
o
l
d
s
t
e
i
nらが用いていた方法であり、体積力点X
b，
t
)を次式で与える
1
7
)
ふ(Xb，t)-V(xb，t)
]
f
(
X
b，
t
}
=αrf~ もかb ， t')-V(Xb ，t')}tt' +β
α
/
，f
J
r
は上述した自由定数であり、それぞれ1fT-， l
ITの次元を有し
)
(
2・2
uはグリッド位
(
X
b，
t
)は境界の速度で、ある.
置での速度 ， V
境界での速度 V
(
X
b，
t
)
=
O，または V
(
X
b，
t
)が定常，式(
2
.
1
)において左辺第一項と右辺第三
項のみ維持すると仮定すると式(
2.
1
)，式(
2
.
2
)より以下の式が得られる.
子
f=α/jqd+伽
(
2
3
)
q=u-V
(
2
4
)
式(
2・3
)は単純減衰振動を表しており ，jが
U
を境界での速度 Yに至らすことを意味し
ている.この方法の欠点は， α
ρ 舟の値が大きいと式 (
2
1
)を不安定にし，時間差分にお
2
2
)
いて非常に細かな時間刻みが必要になる点であるといえる.その改善策として，式 (
の第二項を陰的に扱う式(2~5) が考案されている.
パ
ルt
)叫
(
2
5
)
r
2
.
2 D
i
r
e
c
tf
o
r
c
i
n
g去
体積力の与え方を改良する事により計算の安定性に影響しない方法として d
i
r
e
c
t
t
品
o
r
叩c
i
n
g法が Mo凶
h
叫
1
d
白由定数をうえずず、，計算効率が良く三次元複雑流れへの適応性もある問.これまでにも，
曲線形のノズルから出る渦環構造や球周辺の三次元対称流れ，モーターを備えた IC ピ
ストン内の三次元乱流解析などに応用されている.
d
i
r
e
c
tf
o
r
c
i
n
g法の原理は以下である.式(
2
1
)を陽的に時間差分すると式 (
2
5
)が得られ
る.
un+1_ un
一一一一
d
.
.
t
(
2
5
)
=RHSn +j"
⋮
ここで，右辺第一項 RHSnは式 (
2
.
1
)の粘性項および圧力勾配項から対流項を引いた差で
rにつ
v
川を満たすような Fは式 (
2
5
)を
ある.次ステップにおいて un+1が境界条件 un+1=
n1
γ+1を代入することで単純に
いて解き，境界条件 u
+=
σJ
H
R
7
fJ
(
2・6
)
と与えられる.境界に位置している速度定義点においては， γ叶に境界条件を代入し，
境界内部点においては境界条件を満たすように周辺点より内挿された速度を代入する
'=0とすることで通
ことで体積力 Fが算定できる.境界外に位置している場合は体積力 1
常の N-S方程式を解くこととなる. d
i
r
e
c
tf
o
r
c
i
n
g法の体積力は力学的なプロセスを必要
i
r
e
c
tf
o
r
c
i
n
gといえ，各タイム
とせずに境界で直接的に境界条件を強いるという意味で d
e
e
d
b
a
c
kf
o
r
c
i
n
g法との大きな違いは，/の算定
ステップにおいて境界条件が保たれる. f
の際に自由定数が不要であり，計算の安定性に影響しない点である.本研究では移動境
i
r
e
c
tf
o
r
c
i
n
g法を用いる.
界を扱わないなどの特徴を考慮し，この d
6
2
.
3体積力の内挿方法
速度定義点と境界が一致していれば上述の体積力 fは境界条件をそのまま用いること
が可能だが，一般に複雑境界を扱う場合速度定義点と境界が一致していない.このとき
2
.
2のように
は境界周辺の点の速度から内挿してきた速度を用いなくてはならない.図 色つきの部分が流体領域に設置された物体とし，
0で表されているのは速度定義点であ
る.青Oの速度定義点は流体領域に属しているが，その他の色の速度定義点は物体内部
にある.境界と隣接する点は三通りに分類することができる.まず自分自身が物体内に
属しており上下左右隣接する速度定義点のうち一点のみ流体領域に位置している場合
(図の緑色の点)はその隣接点のうち流体領域に属している速度定義点の速度を用いて
線形内挿が用いられる.隣接する
速度定義点のうち 2点が流体領域
に位置している場合(図の黄色の
O
。
点)は隣接する速度定義点にもう
一点加えた周辺 3点の速度を用い
O
円噛
Q
。
4
て双一次内挿が用いられる.上下
。
O
明
O
左右の隣接点が全て物体内部(図
の白色の点)の場合はその点にお
0
@
l
ー
画
ー
-0
いて境界条件 u=Oを与える.これ
ら内挿された速度を式 (
2
3
) に代
0
0r- -
入することにより物体内部に作
用させる体積力を求める.まず 2
次元形状のときの内挿方法につ
いて述べる事にする.
図2
.
2 内挿の概略図
0
(1)線形内挿
この内挿方法は境界内第一点目の速度と，境界外第一点目の速度を内挿し，境界条件
2
.
3において，色つきの
を満たすように境界内第一点目の速度を求める方法である.図 -
LU
J
今
hL
べ
，
.
，
・
4
一一一
/
/
(
a
)0
くh
h
，の場合
2く
(
b
)h
く
，h
h
3の場合
2く
図2
.
3 線形内挿の概略図
領域が物体領域であり，黒Oで表したのは速度定義点である.赤Oは境界条件を満たす
点である.境界内第一点目の速度は境界から同じ距離のがを用いて
(
2・4
)
U=-u'
と表すことができる(図 2
.3(
a
))
.O
<
h
h1のときがは
2<
u
U
一方 h
<
h
<
h
1
2
3のときがは
九
h一
めることができる.
U1 と境界条件の線形内挿より求
(
2
5
)
U1 と U2から求まる.
Uー
ー
い3 -hタ1 +(h2 - hl炉2
h
3-hl
(
2
6
)
より単純な方法として図 2
.3(
b
)の赤いラインで示したように hl<h2<h3のときも式
(
2・5
)を用いる方法もあるが，その方法では図からもわかるように h
h1の値が格子の分布
2/
などにより発散的に大きな値となってしまい，発散した U により計算の安定性に影響を
h1=100で発散した実例も見受けられる
きたしてしまう.具体的には h2/
h
h1が大きくなると発散が確認された.そこで，流体領域の
2/
1
9
)
本研究でも
U1，U2の勾配を用いて境界
と内挿を行う点までの距離ん離れた地点の速度がを用いることで発散を抑えることがで
きる.
8
(
2
) 双一次内挿
先にも述べた様に上下左右の隣接点のうち 2点が流体領域に位置している場合は，隣
接 2点に l点加えて周辺 3点の速度を用いて，先の線形内挿ではなく双 一次内挿が用い
られる.
干
U2
e
u
ー
l
惨
@ーー砂
Yl
'
‘
L
Z
e
e一
U
一
3
砂
ー Xl
X2
図2
.4 双一次内挿の概略図
2
.4において U は内挿された速度を
図-
U]，U2，U3は流体領域の速度であり，境界
条件は物体の角((図中の赤Oの点)で満たされることとする.このとき物体内の速度 Uは
]，Y2らを用いて
それぞれの速度定義点と境界条件を満たす角までの距離 X]，X2，y
u=一 [(1-αXl-β，
炉 +α(1-β九 +(1-α)βI
J
/(
1
-as)
X，
α =一
一-一一
x，
，+X
'
'2
(
2
7
)
(
2
8
)
)
(
29
由
と表すことができる.挿 一次内挿のときも X] と X2，Y]とY2の比率により線形内挿のとき
と同様に流体領域の 一点目・二点目の勾配を用いることにより発散を抑えることが可能
であると考えられる.また ，Y方向
を用いることができる.
z方向の速度成分 V，W についても同様の内挿方法
(
3
)三次元場での内挿
前述の二次元形状の内挿の場合は X，Yの 2方向のみ考えればよいが，三次元形状にな
ると z方向(流れ直角方向)についても考えなくてはいけなくなる.しかし，内挿が必要
となるのは境界をはさむ面または辺ということになるので二次元形状の内挿方法が応
用できると考えた.三次元形状の場合 X，y，zの各方向の隣接点計 6点が固体領域に属
しており完全に境界内部と認識できる部分を除くと 3種類の境界に分けることができる.
2
.5にその場合わけを示す.色つきの部分が物体領域であり，丸をつけた点は内挿に
図関連する点である.赤丸の点はそれぞれ境界条件を満たす点である.
2
.5(
a
)は物体領域の辺に位置している場合の図であるが，この場合は双一次内挿
図を用いることができる.図 2
.5(
b
)の場合は面に位置している状態であり，線形内挿が
用いられる状態である.そして，二次元形状には見られなかった位置関係が図 2
.5(
c
)
である.この場合は 2次元形状の角の扱いを応用して 3次元双一次内挿の考えを用いる
と式(
2
1
0
)のように表される.
X
l
X
l-sX
lー
か
U=[
(
1
-a -β)
i
Ul+α(1-s
)
仰 2 +(
1
α
)
βi
U3+(
1-a
+
α(
1
-sX1-rル
5 +(
1
α
)
β(
1
-r丸 +αβ (
1
-rル
7]
j(
I
-α
βr
)
4
(
2
1
0
)
、
，
一
X
叶
+
γ一
-一
X
α
(
2
1
1
)
e
(
2
・1
2
)
β =-Y l
Y
I+Y2
一吋
一-
一Z
一Z
'h
一+
γ
'
(
2
1
3
)
4
nHU
l
I
J
I
U2
'
1
42
U1
U2
Us
，
U
y
Y
2
，
X
X2
(
c
)
図2
.5 三次元での内挿の概略図
1
1
第3章
埋込境界法の一般曲線座標系への適応
この章では，実際に本研究で用いた計算プログラムの概略を説明する.運動方程式と
連続式を基本変数である速度のデカルト成分
，圧力 p を用いて解くわけであるが，変
Uj
数の配置は格子の中央に配置されるコロケート格子を用いており，時間進行法は
F
r
a
c
t
i
o
n
a
lS
t
e
pMethod(FS法)を用いたプログラムである.
3
.
1 一般曲線座標での基礎方程式
曲線座標上の点でデカルト座標成分 Ujを未知量とする数値計算では，デカルト座標で
書かれた LES方程式を，境界適合座標に変換して差分計算を行う.デカルト座標 Xjから，
境界適合座標ムへの変換を考慮したとき，基礎方程式は以下のように変換される.
l
δ(JUk) =0
J θC
k
t8
ui
δ￨
(
3
.
1
)
t
8
c
mθC
n8
u
1
ゐ
1
J
t
j
一+一-4Umut+AZV-16mnー」ー -vJ-一一一一一ート
8
t 8
ごm1
勾n
axfaxtθC
ρ
nJ
(
3
.
2
)
運動方程式 (
3
.
2
)には左辺に d
i
r
e
c
tf
o
r
c
i
n
g法により求めた体積力を導入しである.
J
t=detl 竺~l
(
3
.
3
)
lδCm
m
A，
=J- 8x
i
(
3.
4
)
(θCfj
1
u
山
品t
v一
二
J
U
cmn
~~m
1θCmθC
n
J一
一
一
ー
一
一
1
8
x
;8
x
;
(
3
.
5
)
(
3
.
6
)
J
tはデカルト座標から一般曲線座標への変換のためのマトリックスであるヤコビアン
の逆数であり，物理空間でのセルの体積も表している •
Um は速度の反変成分を表してい
る.A
/
ぺ c は座標変換により生じたものである.式(3.1)，式(3.2)は強保存形の保存式
から得られるもので，物理量のセルでの積分値の収支形で表される形となる.
mn
3
.
2 シグマ座標系における微分係数の取り扱い
本研究において扱う座標系は，鉛直方向に直線の座標軸を持つシグマ座標系とする.
コロケート格子を用いた計算では，格子線が互いに直交に交わるものが計算の精度を高
めるとされているが 20) 本研究で扱うような複雑境界を有するの流れ場においては，地
形自体は崖のような勾配が急変する箇所や鉛直方向に近いような傾斜を持ったりしな
いものと仮定し，構造物も一般的な形状である箱型であると仮定できると考える. シ
1
2
グマ座標系を用いる場合，
。
t
x
l
二
BXl=aX3=CX3-0
8C2θC3θご
(
3
.
7
)
8
C
2
という関係式が得られ，数式を部分的に簡略化することができ， 一般座標系と比べ，計
算負荷量を軽減することができるという利点もある.
3
.
3 空間差分法
空間差分において
レイノルズ数が高くなると対流項の扱いが重要となる.なぜなら
ば
， N-S運動方程式の非線形項が高波数成分を次々に生成し，レイノルズ数が高いと，
この高波数成分が散逸されにくくなる.つまり，実用的な格子幅では高波数成分が捉え
きれなくなり，計算格子より短い波長成分が見かけ上長波長のものとして捉えられてし
まい計算を不安定にし，発散に導いてしまうということである.その対策として，本研
究では非線形移流項には梶島の方法 21)を適用し，その他の項には三次精度中心差分を用
いた.
差分計算では，時空間の解像度や差分近似精度次数だけでなく，数値解法，安定のた
めに付加される数値粘性，さらに基礎方程式の「型」によっても計算結果が異なること
が知られている.最近，高精度な数値シミュレーションの要請から
これらに対する検
討が盛んに行われている.その際に，例えば「型」に適合しない差分法において解像度，
精度次数または数値粘性の影響を論じても無意味である. したがって，系統的な信頼性
評価，あるいは高次精度差分の導入の前提として，対流項の型に適応する差分法の選定
が不可欠である.
梶島により，差分の高次化への発展が試みられており，勾配型の対流項について，整
合性の良い差分法，さらにその高次精度化が検討されている.
その検討とは，式(
3
.
8
)のように各速度点において対流項を直接差分化するものである.
ぷ
ぐ十~[:t.~ +v
u
(
3
.
8
)
これは、通常格子の対流型差分と同形式の差分をスタガード格子に適用じたものである.
二次元的な保管が必要な対流速度には一砂を付している.等間隔格子，二次精度では，
+v
1
t
÷
ti-+-fj+vJ+lri+vぺ什
ぺ
u -u
1
j砂 1J=u
[
u3ujdx+vθu
2
'
} 寸
十
V
4
2
h
x
・
U
づ
け
x-
-u
(
3
.
9
)
2
h
となる.勾配型に対する最も直感的な離散化であるが，多次元的な補間を要する点，速
2，.・・)のように広がる点で，スタガード格子にはなじまな
度勾配の差分点は +k
(
k
=1
，
い性質を持つ.この差分法は高次風上を取り入れた乱流の直接計算でも使用されている.
1
3
一方，対流速度の方向に半格子ずらした点で対流項を差分化し計算し，その結果を速
度点に補間するというものも示されている.
(
3
.
1
0
)
トゐ/白 +
v
θu
/砂l
H
j
f
=
l
r
A
u
+
m
v
u
l
J
+
:
J
対流速度の補間および対流項そのものの補聞は一元的に行われ，補間も差分も士
(
k
-1
I
2
)
(
k=1
，
2，・・・)と半格子前後の近接点から展開される. したがって，前述の手法に
比べてスタガード格子との相性が良くなる.この手法は補間法と呼ばれている.
補間法のより本質的な利点は，ゐ/みとか/砂を圧力点 (
i
J
)で連続の式と同じ差分式を用
いて評価することによってもたらされる.つまり，対流項と連続の式との整合性が向上
し，発散型との互換性が高められる.高次中心差分を用いると式(
3
.
9
)と式(
3
.
1
0
)の差は局
所的には速度の発散 (divu) では整理できない.しかし，一般的に差分近似式は格子点
上の値を使った補間式(多項式)の微係数と解釈できるから，連続の式と運動方程式に
おける速度勾配に共通の多項式を用いれば，補間法の原理による勾配型高次差分は発散
型との互換性が増し，高い保存性が期待される.
スタガード格子においては，質量保存が満たされている限り，発散型，勾配型は互換
であり，運動量と 2乗量を保存する性質を持つことがわかる.ただし，この場合の勾配
型の差分化は前述の直接法ではなく補間法でなければならない.
3
.4時間進行法
本研究では Nakayama らの方法
2
2
)を用いて粘性項の対角成分には陰的スキームの
C
r
a
n
k
N
i
c
o
l
s
o
n法を，対流項と粘性項の非対角成分には陽的スキームの二次精度
r
Adams-B
a
s
h
f
o
r
t
h，去を用いる.
(
1
) Crank-Nicolson法
この手法は E
u
l
e
rの陽的差分の式 (
3
.
1
1
)と陰的差分の式 (
3
.
1
2
)の和をとる手法である.そ
3
.
1
3
)では，時間項からくる誤差を互いに打ち消し，時間に関して二次
の結果導かれる式 (
精度となり無条件安定となり精度が良くなる.
+U~ ，
U~+1 - U~
_.U~，' - U~ ，
U~，' - 2u~
一ι一一-'-=-U---'工L 一
一」斗 +V---'工L 一
一一一云
2
企X
d
.
t
n
U
;
. .- U
;
.
n
+
l
n+l
，
'
-u
;
'
'
，
_.U
;
'
ーム一一一，ー =-U---'- ナ
d
.
t
d
.
x
+
n
+
1 司
1
U
;
''
，
-LU;. . U;n+γ
~+v
什
l
¥
d
.
x
2企X
uf+l-uf-Uf
u
;
t
i
l
u
t
i ，U:1- U:
(
3
.
11
)
ιE
L
(
3
.
1
2
)
，
-
L
1i
- 一 一2 一
'2
d
.
t
2l 2d
.
x
d
.
x J
l
:
1
1
2
u
f
I
+
u
Zー Ul
l
2
1
4
7
+
u
l
1i
，V (U
1
4
l
d
.
x2
d
.
x2
(
3
.
1
3
)
式(
3
.
2
)に上述の差分を用いると式(
3
.
1
4
)が導かれる.
J
ι云
ム
=
(3同
f
(
c
;+D ， (u~)+ G (
u
;
)
)
}
(
c
;
t+
D
E
(
uψ
正
qJ/O
、‘.， r
J
3
今
、、.， r
3
弓
，FF
，、、‘..
EEF
、
寸
，.‘、、，，.‘、、
m
、
、utt
m'A
Ui
一
F3-t
j t m δ 一vrb
一一一
一
。
円ο
τ= D A t
c
t
g
ht)m二
云ht)m
二
￡
DE
(
3
.
1
7
)
n
DJ
1
8
)
(
3.
(
3
.
1
9
)
，DE
，GEはそれぞれ粘性項(対
ここで， Qは対流項，1
乙は差分における勾配演算子，D/
角成分)，粘性項(非対角成分)，渦動粘性係数による項である.
(
2
) A
d
a
r
n
s
B
a
s
h
f
o
r
t
h法
nl
を求める外挿型
過去の計算結果 un，uFl，unJを用いて次のステップにおける結果 u
十
ト
の公式であり，修正しない前進する方法である.空間項全体を fとすると，三次精度で
は
，
l
斗 (3gf-gf)
u
f
l
=
u
(
3
.
2
0
)
と表すことができる.
r
a
c
t
i
o
n
a
lS
t
e
pMethodにおける埋込境界法の導入
3
.
5 F
F
r
a
c
t
i
o
n
a
lS
本研究では部分段階法 (
t
e
pMethod) を用いて計算を行っている.この方法
は一つの微分方程式をより簡単な複数の方程式に分割し，それぞれの解の重ね合わせの
原理を用いて元の方程式の解を求める.
3.
3.
)とする.
FS法の段階式にも埋込境界法を適用し，式 (
21
22
)
.，(
p
)
M
(
g
(
O
u)
，
+fi
n1
t
(
g
(
p
n
+
lルf
)
u +=u'+i
1
/ρ
u'=un +
n
2
1
)
(
3.
(
3.
)
22
3.2)を変形すると
ここで gは空間項であり〆はが+1/2の'性格を持った中間値である.式(
1
5
(
I
3
1
刈
(
U
;
n)
_
u
;
n
)=判:
(
c
;
n+
DE
+G (
u
;)
J
lI
2
E
-;ゎ山
レ)ヴイ
(
3
.
2
3
)
)
+GE
(
U
;
l)
}
t
-D
l
1
L2ι=吋 :
r
+
1
)
J
l空
(
3
.
2
4
)
となる.ここで Cは対流項，R
iは差分における勾配演算子，D1 ，DE ，GE はそれぞれ
粘性項(対角成分)，粘性項(非対角成分)，渦動粘性係数による項である.また， (
3
.
2
3
)中
/は次のように求められる.
のf
ff=丘
三竺Lρ-2
i
L
(
c
f
+
D
E
(
ィ
)
+G
A
u
;
)
)
2J-I
(
3
.
2
5
)
+
1
4
(
c
t
nl
+
D
E
b
t
;
?
l
)
+
G
E
(
ィl
)
)
1
乙
(
D
[
(
I
ィ+
'
)
+
D
[
(
ィ
)
)
2J
I
つまり，〆は圧力項無しで運動方程式を満たすものであるから，/も圧力項を除いた関
係式で与える事ができると考えられる.
3
.
2
4
)の φは次式により圧力 p と関係している.
式(
企tT"¥
i
(R
;(
o)
i
W)=(J7117)
( T
l
(
3
.
2
6
)
n1
u+ についての式 (
3
.
2
4
)に φが存在するため ，p の代わりに φについて解く事となる.デ
3
.
2
1
)は Umが定義されるセル表面でも成り立
カルト成分 Uiについての訂正ステップの式 (
3
.
2
4
)を変形すると式(
3
.
2
7
)が得られる.
つ.式 (
(~mn
δr+ i
=U
.
*
.+d
.
.
tl
Gm 一 一 一 l
1
U.~.+I
m
川
(
3
.
2
7
)
)
U
*FJ
m
一r
FM
d
也事
Id
。
円
一
u
Lδ~n
(
3
.
2
8
)
φについてのポアソン方程式は式(3.27)の発散をとり次のようになる.
δf δ
r
+iδU;
円
1
mn
(
3
.
2
9
)
一 一 一 一 -
oCmlδ~n )
〆は H
i
r
s
c
hによる近似因数分解
企tδ~m
2
3
)により反復無しで以下のように求まる.まず(
3
.
2
3
)式
J
.
t3の項が十分小さいと見なし，下式のように展開する.ここで Dの下付きの数字は
を!
計算空間内での方向を示す.
1
(
合引一歩 ]
(
1会吋;*=1
(
合引合 ]
(
1十
三
了D
3
)
u
f
(
3
3
0
)
D2
JD2
lρhu
この式を三段階に分けて解く.
(1-2~~IDl}ir
(
1
+主
，u
l
2J- r=
2J- D
1~l
i
(~'
1
c
A
ui
)
+G~州 ))(cfl+DMl)+G(ufl ， ujl))+ 山1 (
Ui
十
三
占(
i+D
(
3
.
3
1
)
n
)
]
n
n
U
¥illノ
D
U一
V
u
r
t
+
/111¥
¥1ilil-Jノ
D
U一
V
Fi
fell--¥
(
3
.
3
2
)
U
D
μ
一
日
¥11111ノ
U
r
t
+
f'll￨¥
u
一
一
D
¥31111/
Uy
一
r
t
/1111
‘
¥
(
3
.
3
3
)
三重対角行列を解く時の方程式は，
α(
k)
x
(
k-1)+b
(
k)
x
(
k
)+c
(
k)
x
(
k+1
)=f
(
k
)
(
3
.
3
4
)
といった形になる.ここでこの式の係数について簡単に述べると， (
3
.
3
0
)式のご方向の k
番目の点における u
t
'
(
k
)について考えるとき，差分して展開した後に得られる式の
u
t
'
(
k
l
)，U
t
'
(
k
)，u
t
'
(
k
+
1
)にかかる係数がそれぞれ α(
k
)，b
(
k
)，c
(
k
)となり， (
3
.
3
0
)式の
右辺全体が f
(
k
)となる.同様の手順によって三重対角行列を繰り返し解き， u
Jを求める.
次に， u
Jを用いて， (
3
.
2
9
)式の併に関するポアソン方程式を解くのだが，ここではセル
中心のイとセル面上の u;の両方が(
3
.
2
9
)式を満たすようにするために，イから補間され
たu
;を用いる.この時，時間 n+lにおける連続の式が満たされていることを考慮して
いる.これを SOR法を用いて解き，得られる砂川により式(
3.
24
)，(
3.
27
)によってそれぞ
れu
f
+
l，urlを求める.境界での速度を Vとすると境界条件を満足させる体積力 fは次
式で与えられる.
メn
竺
己L
ρ
j去(ct
+DE(ィ
)
+GE(
u
;
)
)
+
i
f(ct-1+E(U 1)+ (U;-l))
fW+l)-jア
タ (U 1)+Dl(
u
t)
了
了
D
n
i-
GE
(
3
.
3
5
)
n
1 i+
fを前章で紹介した内挿法を用いて境界，または境界内で作用させる事により境界条件
を満たすことになる.
1
7
3
.
6計算格子
図3
.1 コロケート格子における変数配置
複雑な地形を扱う際一般曲線座標を用いることはよくあり，本研究においても自然な
起伏のある地表面の境界表現には一般曲線座標を用いる.一般曲線座標を用いる際には
変数の配置に注意が必要である.全ての変数を同一点に配置するレギュラー格子は圧力
分布が市松模様状に振動する計算の不安定を引き起こす可能性もあるため一般曲線座標
にも不向きである.速度と圧力を別々に配置するスタガード格子では計算プログラムが
煩雑になり一般曲線座標への応用が困難である.一般曲線座標では，連続式に速度の反
変成分が用いられることも考慮し本計算にはコロケート格子を用いる
は，基本変数である速度のデカルト成分
2
2
)
2
4
)
この格子で
，圧力 p は計算格子の中央に配置される.連
Uj
続式と移流項の差分には速度を反変成分に変換したものを用いらなければならないが保
存則を離散的に満足させるためにはこの方法が有利である.
1
8
第 4章 数 値 計 算 手 法 の 検 証
4
.
1 計算条件
ここでは埋込境界法を用いた解法の検証を層流について行う.検証に用いるモデル地
形の形状は
(
4
.
1
)
1
+
[
引
で定義される図ー4
.1 に示されるような 2 次元孤立峰である •
高さ方向の座標 ，H は山の高さである.山の頂上の位置
X は流れ方向の座標 ，yは
Xtopは計算領域中央の流入面か
α は斜面のパラメータであり1.1の値をとっている.三次元孤立峰
ら 10Hの地点とした .
の場合のように山の左右からの回りこみは発生しないがパラメータを容易に変化でき，
傾斜を 450'こすることで，丘後
方に比較的大きな剥離領域が
生じるので，剥離流計算の精度
を考察できる.これまでにも実
験
2
5
)や計算 2
6
)が行われている.
計算領域は流れ方向，高さ方向，
5H
奥行きの順に 20Hx5Hx5Hであ
る.傾面の最大傾斜は 450であ
るので境界適合格子では座標
図4
.1 計算領域
軸の交差角は最大 450になり非
直交効果が表れる可能性がある.乱流の計算を行う場合，境界条件として，流入流出条
4
.2のように流入条件
件，壁面境界条件の扱いが重要となる.検証計算では上流端で図 は u=Uo=1
.0，ν=れ =0.0，協=wo=O.Oの
すベり(上面)
U
は主流方向の速度成分 ν は鉛直方
向の速度成分 w は奥行き方向の速度
成分である.また，流出面では自由流
一様流
lF
UO=O'
O
v=Vo=O.
O
w=
Wo=O.O
↓
↓
↓
一様流が流入するとしている.ここで
出条件を，側面には周期条件を，下面
と上面にはそれぞれすべりなし，すべ
すベりなし(下面)
り条件を与えている.一様流の初期状
図4
.2 境界条件
1
9
自由流出
一+
ー+
ー+
態から時間刻み dt
=O
.
O
O
lで，時開発展計算を行い定常状態に達した結果を考察する.こ
のときクーラン数は 0
.
0
0
9
2，格子幅を基にした格子レイノルズ数は 30である.なお，
計算法の検証を目的としているので，流入速度と山の高さを基にしたレイノルズ数は
200と非常に低い値である.計算格子は図 4
.3(
a
)，(
b
)を用い，格子数はどちらも流れ方
向・鉛直方向・流れ直角方向の順に 90X50X25となっている.
'0
5}
}(
5
0
)
刈〉
J
X
(
a
)境界適合座標系
5
H(
50
)
(
b
)直交座標系
図4
.3 計算格子
2
0
検証に用いる手法のうち，境界適合座標系は丘地形に沿って作られた計算格子状で計
22
) により， 2次元形状の流れ場での層流，
算が行われる.この手法は， Nakayamae
ta
.
1
乱流に対して検証が行われており，結果を比較することにより計算手法の妥当性を検証
4
.4のように各セル中心点が地形内に属
することが出来ると考えた.ー矩形近似法とは図 するか，流体領域に属するかの判断により地表面の境界表現を行う手法である.この手
法は格子生成が容易であり，物体境界の処理も容易であり，直交格子の利点を生かした
まま計算を行うことが可能である.しかし，地表面は連続的に変化しているが境界条件
は定義された点でのみ与えられるため，精度上に問題が残ると考えられる.矩形近似法
と比べることでも計算手法の妥当性を検証することが出来ると考える.
I I I I
J
一
実際
、
界
r-li-lin
十r
『見
TI
口流体セル
、
-、
、
.
レ
./
、
'
-
レ/
卜¥
./
v
./
一
、
、
¥
、
、
斗¥
~い/
LJ
物体セル
、
句、
、
f
IIII
i
I
I
I
I
I
I
I
I
I
I
I
I
口『
図4
.4 矩形近似法の概略図
4
.
2 検証結果と考察
4
.5(
a
)は山の頂上から 見て 2H風上側の地点での主流方向速度の鉛直分布の結果で
図ある.実線は境界適合格子を
x印は矩形近似法を， 0印は埋込境界法を用いた計算結
果を x/H=8.0 の鉛直面に内挿したものである • x
/H=8.0は山の風上で減速している位置で
4
.5(
a
)から見られるように全ての結果は良好に 一致している.図 4
.5(
b
)は x/H=lO
図の位置であり，丘の頂上での流速分布である.頂上において加速していることが確認で
4
.5(
c
)(
d
)は山の頂上から見て H，
きる.ここでも全て計算結果はよく一致している.図 2H 風下側の地点での主流方向速度の鉛直分布の結果である • x/H=11
.0は山の風下で逆流
4
.5(
c
)から見れるように境界適合格子と埋込境界法による結果は逆
している位置で図 -
4
.5(
d
)は
流が確認できるのに対して矩形近似法では逆流が確認できなかった.図 x/H=12.0での結果である.全て逆流が確認できるが，y/H
亘2において矩形近似法は流速
が大きく出ていることが確認できる.これらの差異も風土側のときと同様に丘地形が階
段状になってしまうことが要因と考えられる.矩形近似法により実際の丘が階段状にな
り，実際の高さと最大で O.lHの違いが確認された.
2
1
5.
00
5.
00
一一一 境界適合格子
× 矩形近似法
o IBM
4
.
0
0
4
.
0
0
3.
00
3
.
0
0
と
足
、
、
ロ
コ
為、
戸ヘ
¥、
2.
00
2.
00
1
.00
1
.00
0
.
0
0
0
.
0
0
0
.
0
0
吊
0
.
0
0
1
.5
0
1
.0
0
0.
50
u
/
Uo
1
.5
0
(
b
)x
/
H=1
O.
0
(
a
)x
/H
=8
.0
5
.
0
0
1
.00
0
.
5
0
u
/
U
(
)
5.
00
グ
ヨ
4.
00
4.
00
3.
00
3.
00
:
t
¥
、
巴
と
¥
h
岸、
2.
00
2.
00
1
.00
1
.
00
0
.
0
0
-0.
50 0.
00
0.
00
-0.
50 0.
00
0.
50 1
.
00
u
/
Uo
1
.
50
0.
50 1
.
00
u
/
U
(
)
(
d
lx/
H=12
.
0
(
c
)x/
H=11
.0
図4
.
5 流速分布図
2
2
1
.
50
5
1
.4
0
.74
ylH
O L-
ニLー
て
;
;
?
'
L_~
(一
え‘
、ー
O
「
1
0
20
xl
H
(
a
)境界適合格子
5
ノ
卜￨
。~
/
卜
￨
ヘ
/
/
ylH
，
OL
ー
<;-一一
==エニニ-
l
一
ー一
「
1
0
O
20
xl
H
(
b
)矩形近似法
5r- I
~ . 74
yl
H
I
1
OL- I
~
-
，
= ，ー
るー
O
考、
さ
、き
、ー
ー
】
1
0
、
.
、
‘
2
0
xl
H
(
c
)I
B
M (埋込境界法)
図4
.6 流速等値線図
4
.6は流れ直角方向中央断面伝I
H
=
2
.
5
)における流速等値線図である.この結果も先
図の流速分布と同様に時間発展計算を行い定常状態に達した時の結果である.境界適合格
子を用いた結果と埋込境界法を用いた結果に比べて，矩形近似法を用いた結果では丘後
方の剥離領域において差が生じていることが分かる.矩形近似法では丘の高さに実際と
の差が生じ，結果として丘後方の逆流領域が小さくなったのではないかと考えられる.
境界適合格子を用いた結果と埋込境界法を用 いた結果は良好な 一致が見られることも
わかる.
2
3
5
・・
-・
・
・
・・
・ 一一
・
・0.75
F
~. 1 7
y/H
lm
~ケプアミ二
二
O
ー
ー可
1
0
O
x/H
(
a
)境界適合格子
5
0
.7
5
y/H
O
1
0
O
20
x必7
(
b
)矩形近似法
5
y必7
O
O
1
0
20
x/H
(
c
)I
B
M (埋込境界法)
4
.7
図-
z方向渦度等値線図
図4
.7は Z方向渦度等値線図である.丘の頂上付近において境界適合格子，埋込境界
法の結果と比べて矩形近似は値が若干大きくなっている.これは地形が階段状になり，
直立する壁に当たることが原因ではないかと考えられる.
以上の検証より，矩形近似法を用いた結果と境界適合格子を用いた結果の比較では矩
形近似法において差異が見られたのに対して，埋込境界法を用いた結果と境界適合格子
を用いた結果の比較では結果の良好な一致が確認でき，本研究の計算手法の妥当性が証
明されたといえる.
4
.
6，図 4
.7において埋込境界法の結果の x/H=8.0付近で計算の不安定が確認でき
図2
4
る.これは差分法に中心差分を用いたり，不等間隔格子において格子幅が細かくなる位
置で確認できるものと考えられる.後の 5章に示す結果において流れ方向の格子を等間
隔にした際この不安定が無くなったことから，格子幅の影響であったのではないかと考
えられる.
2
5
第 5章
一般座標系 IBMによる 2次元障害物のある地形上気流解析
第 4 章において埋込境界法を導入した計算法の妥当性は示された.そこで本章では，
そのプログラムを用いてモデル地形上に存在する構造物を考慮した流れ計算の結果を
示す.
5
.
1 計算対象
地形には余弦の 二乗で表すことのできる二次元丘地形を用いた.余弦の二乗で表され
る地形は，解析に用いる代表的なモデル地形のひとつ
2
6
)でもあり，検証で用いた丘地形
よりもよりなだらかな斜面を確保でき，実際の丘に近い形状であると考えた.主流方向
を X，鉛直方向を y とする時，地形は次式 (
5
.
1
)と表される.
吋
Y
C=
2
π
王
子
)
7
.
0三
五x壬1
3
.
0
(
5
.
1
)
ここでんは地形の標 高，H は丘の高さ， λは底面の長さで λ
=6Hである .形状の断面は
.1になる.山の頂上の位置 Xtopは計算領域中央の流入面からlOHの地点とした.
図
ー5
y
一
一
一
一
一
一
'
ー
Uo
一
一
一
一
ー
一
一
一
一
一
一
'
ー
x
6H
図5
.1 計算に用いた
2次元丘地形
5
.2に示した奥行き
地形上に設置する構造物は図 -
y
無限長の直方体を用いた.ん，んはそれぞれ主流方向，
鉛直方向の構造物の長さである.
このん，んを任意に変化させることにより構造物の大
~ð ~~{tðit
矩…のーと比附
う.計算格子は図5
.
3のような格子数が流れ方向・
l
鉛直方向・流れ直角方向の順に 90X50X25となって
z
いる.一様流入の速さ Uoと山の高さ H で無次元化し
図5
.2 2次元構造物
2
6
打
H
た時のレイノルズ数 Re=~ O~~
は 100 としている.
V
51
吋;
J
i
J
l
図5
.3 計算格子
UnH
レイノルズ数 Re=~ O~~ =100の時，地表面近傍の鉛直方向格子間隔は O.IH以下である
V
ことが要求される
2
5
)
本研究では鉛直方向に等比数列的な不等間隔格子を用いており，
地表面近傍において鉛直方向の格子幅は 0.02Hである.境界条件は計算方法の検証の時
と同様に様流入，自由流出，地表面すべりなし，左右両端面と上端面はすべりとし
ている.
5
.
2計算結果と考察
境界と速度定義点が一致している場合は矩形近似法でも速度定義点において境界条
件を与えることができるため，正確に境界を表現することが可能である.そこで，まず
はじめに境界と速度定義点が一致している場合の計算結果を示すこととする.一致して
いる時は矩形近似法と埋込境界法を用いた計算結果は一致するはずである.その後，任
意の位置に境界がずれた時に両手法の計算結果にどのような差異が生じるかを示して
いくこととする.
2
7
1
) 境界と速度定義点が一致している場合 (
c
a
s
e
l
)
.
4に示すとおりである.拡大図において格
c
a
s
e
l における地形と構造物の配置は図ー 5
子の交点が速度定義点としてある.計算格子の形状の都合により矩形近似法と埋込境界
法において y方向の境界(構造物の天井の高さ)は速度定義点と完全に一致させることは
できないがほぼ一致するように設置しである. x方向の境界は速度定義点と一致させて
ある.
x
7
8
9
10
xl
H
図5
.4 構造物配置図
1
1
1
2
13
(
c
a
s
e1
)
図5
.5は z方向中央断面における流速の等値線図である.両結果を比較しても丘と構
造物を越えるとき加速し，構造物後方において逆流領域が形成されていることも見て取
れる.この結果においても矩形近似法と埋込境界法の結果はよく一致していることが分
5
.6は z方向渦度の等値線図である.両計算において構造物の上部において絶
かる.図 対値の大きい渦度が確認できる.以上の結果より，境界と速度定義点が一致している場
合は矩形近似法と埋込境界法において同様の結果が得られる.
4章での検証において埋込境界法の結果に x
1H=8付近において数値振動が確認できた.
5
.
5，図 5
.6を見ると，その振動がおさまっていることが確認できる.図 5
.
3
しかし図 に示した通り解析に用いた計算格子は流れ方向に等間隔となっているため，格子間隔が
狭まる地点が無くなったため，振動もおさまった結果が示されたと考えられる.
2
8
5
yメH
O
(
a
)矩形近似法
5
yメH
O
20
x/
H
(
b
)I
B
M (埋込境界法)
図5
.5 流速等値線図
(
c
a
s
e
l
)
•
5
0
.82
←
一一
一
一
ー"ー
一 ・ ー
ー
ー
吾与 主ε
←
-一-~一
，
ー
、
A一
、
、一
ー
町
一
ー
‘
白、
、
、
、
一
1
・
0'-----
、
ー
ー
、
一
一
A
>.
7
，
.
、一
一
、
、
‘
ー
、
一
、
、
r
‘
:
;'
-ー'
e
10
20
xlH
(
a
)矩形近似法
•
5
0.82
川ト
6.
7
I rLI~
/
、
.
.
:~，
、
・一. ~
0'------
10
20
x/
H
(
b
)I
B
M (埋込境界法)
図5
.6 z方向渦度等値線図
2
9
(
c
a
s
e1
)
2
) 境界と速度定義点が一致していない場合 1
(
c
a
s
e
2
)
Case2における地形と構造物の配置は図ー 5
.
7に示すとおりである.構造物の大きさは
ん=
l
.OH
，ん
=O.5Hであり ，x=12.0Hに構造物の前面がくる様に設置した.
y
EEEa
l mg
(
c
a
s
e
2
)
x
/
H=1
2
-入
今
J
EEA
z
図5
.7 構造物配置図
x
x
/
H=1
3
1
.00
100
0.
80
080
，
×
，
0.
60 ・
¥ζ
ー
、
040
，
:
>
.
060
三
040
，
，
020
020
000
-0，010
0，00
-0，010
，
，
，
-0，005
0，000
0，005
0，010
-0.
005
(
a
)境界前面
0，005
0.
010
(
b
)境界背面
5
.8 境界での流速分布
図図5
.
8 を見ると，
0.
000
u
/
Uo
u
/
Uo
(
c
a
s
e
2
)
(
a
)においては速度定義点と境界位置が一致していたため両計算結
果に差異はない.しかし
(
b
)において埋込境界法の結果は境界条件
u
=
Oを満たしてい
るのに対して矩形近似法では負の流速が確認できる.これは矩形近似法において境界条
件を満たす点が -x方向にずれてしまい，構造物後方の逆流領域の流速をプロットしてい
(
i
，
j
.
k
)と x
(
i
+
lふk
)の値をとるとき，境
ると考えられる.矩形近似法では境界の x座標が x
(
i
，
j
.
k
)を認識するためこのような差が生じるのである.図 5
.9は構
界の x座標としては x
造物周辺における流速ベクトル図に速度 0の等値線図を加えたものである.図中の赤い
領域が構造物であり，青い線は速度 0の等値線を示している.この図からも前面の壁で
は両結果とも境界条件を満たしており，その一つ手前の速度定義点では速度が発生して
いるのがわかる.しかし，後方においては矩形近似法において境界位置がずれているこ
とが分かる.図からわかるように矩形近似法では速度定義点を結ぶ形で境界が表現され
ている .y方向の境界においても埋込境界法の結果は内挿により y
/
H
=
O
.
5において速度 0
3
0
の等値線が見られるのに対して矩形近似法では y/H
孟0
.
5における最も大きな y座標を持
つ点を連ねられており ，正確に境界高さが表現されていないことが分かる.
一一~
~
J
ー-ー 官房
¥、
ーー
一 一一.
二
，
、
、、
、
¥
骨一
、
、、h
、、ー-
」
0.
5
¥
、
¥
事 ¥
『ト』、
y
l
H
司 、『
司
ト
『
ー
司トー・
。
‘
・
ー
ー
1
2
1
3
x
l
H
(
a
)矩形近似法
一一~
一
一
骨
e
.
.
、
ご、
ー『
ー
-一
一
一、
.
J
司 『、
a戸
、¥-
司 、、
、¥
、
0.
5
--
/
1
y
l
H
。
1
2
x
l
H
1
3
(
b
)I
B
M (埋込境界法)
5
.9 構造物周辺の流速ベクトル図と速度 Oの等値線図 (
図c
a
s
e
2
)
31
0.
5
y
l
H
。
1
2
1
3
x
/
H
(
a
)矩形近似法
一
一
一
一
一一ー
ーー
』ー
一
一
ー一ー』ー
ー一ー一、ーー一ー一一一一ー-~
0
.
5
y
l
H
。
1
2
1
3
x
l
H
(
b
)I
B
M (埋込境界法)
図5
.1
0 構造物周辺における流線図 (
c
a
s
e
2
)
3
2
図5.1
0は構造物周辺の流線図を同じ条件で描いたものを表したものである.矩形近
似法において構造物が実寸より小さく表現された結果，構造物後方で流線が入り込んで
しまっていることが確認できる.
3
) 境界と速度定義点が一致していない場合 2
(
c
a
s
e
3
)
C
a
s
e
3における地形と構造物の配置は図 5
.1
1に示すとおりである .構造物の大きさは
ん=2.0H，l
y
=O.
75Hであり， x=1
1
.5Hに構造物の前面がくる様に設置した.
y
x
図5
.1
1 構造物配置図 (
c
a
s
e
3
l
x
/
H=1
1
.5
x
/
H=1
3
.
5
1
.
00
×
0.
80
0
.
8
0
d
>X
.'
"0.
60
_"
，
三 0.
40
×
〈
0.
20
0.
00
-0
.
0
1
1
.
0
0
-0.
005
I ~ 矩形近似法
入
oI
B
M
0
0.
005
0.
60
0.
40
0.
20
0.
0
1
0.
00
-0.
010
-0.
005
u/
Uo
0.
000
0.
005
0.
010
u
/
Uo
(
a
)境界前面
(
b
)境界背面
5
.1
2 境界での流速分布 (
c
a
s
e
3
)
図図5
.1
2から埋込境界法の結果は両境界において境界条件
u
=
Oを満たして いるのに対
して矩形近似法では流速が確認できる，これは矩形近似法において境界条件を満たす点
x方向に，背面では x方向にずれてしまいこのような差異が生 じるからであ
が前面では +
5
.1
3は構造物周辺における流速ベク トル図に速度 0の等値線図を加えたもので
る.図 ある.この図より矩形近似法において境界条件を満たす点が前面，背面においてそれぞ、
れ構造物を実寸より縮小する方向にずれている.埋込境界法による計算結果は結果は y
方向の境界においても内挿を行い，境界条件がよく満たされていることが確認できる .
5
.1
4は構造物周辺の流線を同じ条件で描いたものであるが，こち らにおいても構造
図物の 上方において，矩形近似法の結果は構造物に入り込むことが確認できるが，埋込境
3
3
界法の結果は構造物の高さを正確に認識できており入り込みは見られない.
1r
0
.
7
5
-
ー~-一一-
~
d
d
，
、
戸-IA
、
崎
EE-
・
・
。
x/H
(
a
)矩形近似法
0
.
7
5
司
"
-
、
守、民
J
.
，
".r
Y/H
.
.
-"
。
.
.
.
炉
』
1
3
.
5
1
1
.5
x/H
(
b
)I
B
M (埋込境界法)
c
a
s
e
3
)
図5
.1
3 構造物周辺の流速ベクトル図 (
3
4
-JFV4dlh
k¥¥¥
ーー噌~
.、司司﹃司
¥¥¥¥と﹄
y
/
H
一
一
ー
一
、
一
一
『
一
一
一
一
、
ー
ー
0.
75
y
l
H
O
11
.5
x
l
H
1
3
.
5
(
a
)矩形近似法
0.
75
y
l
H
O
1
1
.5
x
l
H
13.
5
(
b
)I
B
M (埋込境界法)
図5
.1
4 構造物周辺における流線図 (
c
a
s
e
3
)
以上の結果より，境界と速度定義点がずれた場合に矩形近似法では正確な境界を表現
することは不可能である.直交格子を用いて，境界位置に格子が来るように格子形成す
る方法も考えられるが，直交格子を用いることにより地形を解像することが不可能にな
ってしまう.しかし埋込境界法を用いることで，任意の位置の境界を比較的正確に表現
することが可能であることが分かった.
3
5
第 6章 一 般 座 標 系 IBMによる 3次元構造物のある地形上気流解析
5章では地形上に 2次元の構造物を設置して解析を行いその結果を示してきたが，実
際に存在する構造物は 3次元形状であり，流れ直角方向においてもそのサイズは有限で
ある.その際用いられるのは 3次元内挿法である.この章では 3次元構造物を地形上に
設置し，解析した結果を示していく.
6
.
1
計算対象
6
.1に示し
用いる地形と計算格子は 5章と同様である.地形上に設置する構造物は図 zはそれぞれ主流方向，鉛直方向，流れ直角
た奥行き有限長の直方体を用いた.ん，ん.l
方向の構造物の長さである.この lx. ん • l
zを任意に変化させることにより構造物の大き
さを変化させ，矩形近似法の計算結果と比較を行う.
一様流入の速さ Uoと山の高さ H で無次元化した時
y
打 H
のレイノルズ数 Re=~ O'ーは 100 としている. 境
V
界条件は同様に， 一様流入，自由流出，地表面すべ
りなし，左右両端面と上端面はすべりとしている.
6
.
2
計算結果と考察
ら
1
)境界と速度定義点が 一致している場合 (
c
a
s
e1
)
x
図6
.1 3次元構造物
丘の頂上に境界と速度定義点が一致するように
l
zを設定し構造物を設置した(図 6
.2
)
. こ
ん×ん X
y
3.0H
~)
1
と
O.
67H
O.67H
z
(
a
)x-y平面図
(
b
)x-z平面図
6
.2 構造物の配置図
図-
(
c
a
s
e
l
)
のケースでは矩形近似法と IBM の結果は良好な一致を見せており構造物に入り込む流
6
.
3
).後方に構造物後方で発生した剥離渦も両計
線はどちらの結果も見られない(図 3
6
算結果ともに良好に一致している.図 6
.
4，図 6
.
5はそれぞれ構造物周辺における y方
向及び z方向の渦度等値線図である.図 6
.4をみると構造物手前の左右において渦度の
6
.5でも両計算結果は良好に一致しているこ
高い部分が両結果において確認できる.図 とが確認できる.
弓
F
ー
一戸一-
ZIH
『
・
ー
ー
目、
-ー
ー
、
ー
。
U
1
0
XI
何
(
a
)矩形近似法
z
戸
H
O
1
0
x
ぷf
(
b
)I
B
M (埋込境界法)
6
.3 構造物周辺における流線図
図-
3
7
(
c
a
s
e
l
l
5
0.
3
4
0.
34
.
、
.
gチ
f
ノ
。
一一6
1
6
1
0
x
l
H
(
a
)矩形近似法
5
。f
‘ー
0.
3
4
0
.
34
Z
I
チf
、-
。
6
1
0
xl
H
(
b
)I
B
M (埋込境界法)
図6
.4 構造物周辺における
3
8
y方向渦度等値線図 (
c
a
s
el
)
1
6
0
.
5
2
y
グ
子
、
ー
O
1
6
1
0
x
l
H
(
a
)矩形近似法
0
.
5
2
y
l
何
唱
-11 4U
EEa-
O
1
0
x
l
H
(
b
)I
B
M (埋込境界法)
6
.
5 構造物周辺における Z方向渦度等値線図
図-
3
9
(
c
a
s
e
l
)
c
a
s
e
2
)
2
) 境界と速度定義点がずれている場合 (
丘の後方の剥離領域に構造物を図 6
.6 のように設置した • I
x
XんXI
z
=1
.0HX0
.
5H X1
.0
H と設定し，境界と速度定義点は一致していない.
﹂
-
y
1
.0H
1
.0H
z
1
.0H
(
a
)x
y平面図
(
b
)x
z平面図
図6
.
6 構造物の配置図 (
c
a
s
e
2
)
図6
.7は地表面から約 0.2Hの高さの j=4の面における構造物周辺の流線図である.矩
形近似法の図において，構造物の領域を示す紫の領域内に書かれた四角は矩形近似法に
おいて認識された構造物の形状である.構造物が逆流領域に設置されているため図の右
から左への流れである.どちらの計算結果も構造物に対して左右対称な流れが形成され
ているが，矩形近似法では設定した境界位置をとらえることができずに構造物のサイズ
が小さく表現されてしまい流線が構造物内に入り込んでしまっていることが確認でき
る 一方，埋込境界法の結果においては速度定義点聞に存在する任意位置の境界を表現
できており，実寸の構造物をよける流れが確認でき，構造物内に流線が入り込むことは
6
.8は地表面から約 O.
2
H の高さの j
=
4における構造物周辺にお
ないことが分かる.図 ける y方向渦度等値線図である.こちらの結果も構造物に対して左右対称な流れが形成
されているが，構造物の大きさに差異があるため渦度の発生にも違いが見られる.
埋込境界法において内挿を用いて任意の大きさの構造物を表現する際に，構造物内に
どの程度の速度定義点があればよいか
1 点のみ存在するような状況では正確に構造物
の幅を認識することはできない.それは，内挿によってその一点のみが用いられるため
である.最低 x，y，z各方向に 2点ずつ格子点が存在すればその大きさをとらえる事が
可能と考えられる.
4
0
5
ZIH
ー'
ーー
-
-
l
O
o
o
-
-EA
l
XIH
(
a
)矩形近似法
5
ZIH
-EA
A
00
l
'EE
﹄
l
O
XIH
(
b
)I
B
M (埋込境界法)
図6
.7 構造物周辺における流線図(j
=
4の面) (
c
a
s
e
2
)
41
5
0
.3
5
・0
.
3
5
ZIH
O
x
1
H
(
a
)矩形近似法
5
0
.35
0
.3
5
gチ
子
O
ルデT
(
b
)I
B
M (埋込境界法)
図6
.8 構造物周辺における y方向渦度等値線図(j
=
4の面) (
c
a
s
e
2
)
4
2
3
) 構造物が複数存在する場合 (
c
a
s
e
3
)
丘の後方の剥離領域に構造物を図 6
.9のように 2個設置した.構造物の大きさはそれ
z
=1
.0HX0
.
5HX1
.0H と設定した.
ぞれん×ん X/
1
.0
1
.0
己
1
.0H
x
(
a
)鳥撒図
(
b
)x-z平面図
図6
.9 構造物の配置図
(
c
a
s
e
3
)
図6
.1
0は構造物周辺における流線図である.両計算結果ともそれぞれの構造物周り
a
s
e
2と同様に矩形近似法においては構
において左右対称な流れが確認できる.しかし c
6
.1
1 は構造物周辺における y方向渦度等値
造物内に流線の入り込みが確認できる.図 線図である.こちらの結果も構造物に対して左右対称な流れが形成されているが，構造
物の大きさに差異があるため渦度の発生にも違いが見られる.
複数の構造物が存在する場合，その大きさとともにその距離と，その聞に存在する速
度定義点の数に注目すると，異なる構造物の聞には最低一つの速度定義点が必要となる
と考えられる.
4
3
5
.~ーー一ー
z
1
H
l
O
z
唱
--a
-A
-
1
7
xJ
H
(
a
)矩形近似法
5
ー
、ー
ー』
ー
ー
ー
ー
ー
ー一一
一
z
1
H
- 一一一一一ーー
一
ーー一
-
nU
A
唱E 4
・
・
・
唱
L
がH
(
b
)I
B
M (埋込境界法)
図6
.1
0 構造物周辺における流線図 (
c
a
s
e
3
)
4
4
1
7
/lilih
〆 /{¥
5
zlH
111
O
A
-EE
I
司
x/
H
(
a
)矩形近似法
5
0.
1
zlH
。
がH
(
b
)I
B
M (埋込境界法)
図6
.1
1 構造物周辺における y方向渦度等値線図 (
c
a
s
e
3
)
4
5
0
.
1
内挿の際に，構造物の角はとがっており，特殊な形状と認識できるが，検証の結果を
見てみると良く表現されていることが分かる.以上の結果より，境界適合格子に埋込境
界法の 3次元内挿を施すことで，矩形近似法では不可能である境界適合格子上での 3次
元物体の境界を表現する手法を開発できたと考えられる.
4
6
第 7章.一般座標系 IBMの実地形上気流解析への適用
5，6章では，本研究で作成された一般座標計 IBMのプログラムの妥当性を，モデ
ル丘地形に構造物を 設置した流れに適用することにより明らかにした.この章では，
この計算手法の実地形への適応を行った.
7
.
1 計算対象
今回用いた実地形は富山県にある山岳地帯である.地形の概略と 計算格子を図一 7
.1
に示す.地形は標高も一緒に示されている.
800
主/
NL.
図-7.1 計算領域と計算格子
地形は国土地理院が発行する数値地図を基に作成
27)しており，データの水平方向の
解像は X，Z方向共に 50m間隔であり ，y方向は 1mである.計算格子数は X，y，z方
1
9である.格子間隔は，各方向共に等比数列的な不等間隔と
向にそれぞれ 117，30，1
している.構造物は地形のスケールを考慮しサイズを仮定した. 計算領域で無次元長
さ 1が実際のスケールでは 350mであり，構造物は 一辺が計算格子状で 0
.
2から 0
.
3，
47
N
て士一
図7
.
2 構造物の位置
.
0
4から 0
.
0
6とした.実際のスケールでは約 100m四方で高さが 30mか
高さ方向に 0
ら 40m.ホテル並みの大きさである.これらの建物を図 7
.2に示した位置に配置した.
建物設置付近はスキー場が実在しており，実際に大型構造物が存在してもおかしくな
い地域である.計算の結果は，建物周辺(図中の枠内)において，建物を設置しない時
と設置した時の流れの違いを見ていくことにする.境界条件はこれまでの検証と同じ
とした.本計算では北風が流入条件として与えられていることとなる.
7
.
2 計算結果と考察
図7
.
3，図一7
.
4，図一7
.5はともに風速ベクトルとスカラー風速コンター図である.
図中において建物は黒い四角で表している.図-7.3は地表面から約 5メートル付近の
高さである.四角で囲ったところが建物でありコンター図と流速から内挿により境界
が表現されて建物を認識できていることが分かる.建物がない場合と比べても建物前
方において減速が見られ，建物の横に回りこみ，建物の聞を通る際に加速が生じてい
ることも確認できる.函-7.4は地表面から約 30メートル付近の高さであり，建物の
.
3の結果より流速が早く，
およそ中心ほどの高さである.全体として地表面付近の図ー 7
建物も良好に表現されていることが分かる.図-7.5は地表面から約 50メートル付近
の高さであり，建物の高さと同じかそれ以上の高さの地点である.図の一番下にあっ
た構造物の位置では建物の上空になっており，建物がない時の結果と一致している.
48
0.
06
0
。
0
N
.~
U
ー
ー- =0.
1
(
a
)構造物がない時
06
入r
ヲ
u
ー
一
ー =0
.
1
(
b
)構造物がある時
7
.3 j
=
2における風速ベクトルとスカラー風速コンター図
図49
λ
r
m
b
~-1. 0
(
a
)構造物がない時
0.
2
0
J
九
f
~~
1
.0
(
b
)構造物がある時
=
4における風速ベクトルとスカラー風速コンター図
図-7.4 j
50
N
寸
“
(
a
)構造物がない時
N
(
b
)構造物がある時
7
.5 j
=
7における風速ベクトルとスカラー風速コンター図
図5
1
2.
4
2
.
2
戸
グT
2
.
0
--一一ー一
一一ー
l
.8
1
.6
6.
0
x
1
H
8.
0
7.
0
図7
.6 構造物を通る中央断面における流線図
2.
4
2.
2
0
.
1
0.
1
y
/
H
2
.
0
l
.8
ー一
ー一
戸
一ー
ー
一ー一
一
一
、
d岡
戸
ー
、
、-
、
1
.6
町、
ー、
6.
0
x
1
H
7.
0
8
.
0
図7
.7 構造物を通る中央断面における z方向渦度コンター図
52
、
(
a
)
構造物がない時
、
気
(
b
)構造物がある時
図
ー7
.8 北東上空から見下ろした地表面付近の流線図
5
3
図7
.
6，図-7.7は図-7.2に示した北から 2番目の構造物の x
y中央断面での流線図
と z 方向渦度コンター図である •
x
y平面においても構造物を認識した流れが確認で
きる.
.8は上空から見下ろした時の流線図である.構造物がない時は地
最後に示した図ー 7
形の起伏の影響のみを受けて流線が曲がったりするのに対して，構造物があるときは
構造物の影響も流線に反映されていることが分かる.
以上の結果より，本計算プログラムは実地形上においても建物の位置を任意に設定
でき，その大きさも自由に変えられることが分かる.しかし，実現象への応用となる
と
， Reが高い乱流解析となる.その場合は地形の境界条件に対数則などのモデルを導
入することでより実現象により近いかたちでとらえることが可能となると考えられる.
54
第8章埋込境界法の一般化
7章までは，地形上に存在する構造物等の境界条件 (
n
o
s
l
i
p条件)を埋込境界法によ
り表現する手法の検証等を行ってきた.しかし，この埋込境界法の可能性を運動方程式
以外の熱や物質拡散の輸送方程式にも拡張できれば，河川や沿岸への温水や混合物質の
拡散の予測法の改良につながると期待される.
この章では，温度や濃度など速度場だけでなく Neumann条件やその他の条件に適用で
きる埋込境界法を提案し，計算例を示すことによりその適用性，実用性などを検討する.
この方法で一般化した埋込境界法は乱流エネルギーや散逸などの乱流統計量の計算にも
応用可能である.
8
.
1Neumann型境界条件への適用
速度ベクトルの粘着境界条件のように値がゼロとなるような場合はこれまで検証し
てきた通常の埋込境界法の適用が可能だが，その他の境界条件の場合は基礎式を導き直
す必要がある.ただし，この場合も境界条件は
un+1=un +M(RHS+f
)
(
8
.
1
)
と外力 fを与えることにより満たすようにする.速度境界条件に境界法線方向勾配を設
定する Nuemann条件
:
:
L
=
v
日
生
(
8
.
2
)
8
n
を課す場合を考える.ここで nは境界外向き垂直方向距離，内はそこで与えられた速度
の法線方向勾配である.この場合 fを
θ f θ RHS
dun
V
'
h一一一)/
!
:
J
.
t
一一=一一一一一 +(
θ n δn
d
n
(
8
.
3
)
x
x
z
u
Ill111
﹂
寸
Aa
u
Ee
¥yiij'j
n
，
d
V
9
J
ffi--¥
F
・
+
、
m
ち﹁卜L
す
わこ
な/
値
J
V
主
ナ
つ
よ
よ
寸/
る
と
(
8.
4
)
8
.
1
)に代入し境界に対する垂直方向微分をとり境界位置 Xbで評価する
とする.これを式 (
と
5
5
θun+1
δn
L
1
{
苧
+
号
)
=
ぞ+
=
千
件
+
{
乎+
b引
(
8
.
5
)
=v'b
となり Neumann条件(式(
8
.
2
)
) を満たすことになる.数値計算では fは離散化したデル
タ関数で近似することになるので格子点 Xjでの値は
/Llt~(Xi
f
(
X
i
)=l-RHS州 +(
V
'
bn
i-Un)
一喝)
(
8
.
6
)
となる.境界での速度勾配 v
¥がゼロの場合は
中川 Xb)
↓
n
f(Xd= RHS(Xi)-u /
(
8
.
7
)
である.壁画で速度を指定するのではなく応力を指定する場合(壁面応力モデルなど)
や粘性応力を指定する場合 fはその応力テンソルの発散になる.
8
.
2運動方程式以外の輸送方程式への埋込境界法の応用
前節では運動方程式より流速についての時間進行解を求める場合を示した.対象が
運動方程式でなく他の輸送方程式の場合もほぼ同様な方法が適用できる.
温度 T
，混合物濃度 C についての輸送方程式に適用する場合，外力に相当するで熱輸
送 gおよび物質輸送 hを付加し
θT
τ
=
-= -¥
1.
(
uT)+β¥
1
.
¥
1T+g
(
8
.
8
)
C
J
t
θC
(
8
.
9
)
τ
=
-= ¥
1
.(uC)+y¥
1
.¥
1C+h
o
t
を解けば良いことになる.ここで s， ァは熱伝導係数，物質拡散係数である.この場合
g，hは
仰 )=(-RHST+(
九
n
T)
叫仰
X
b
)
(
8.
10
)
n/企t
抑)=(一問先 +(C
)仰 -Xb)
b-C )
(
8
.
1
1)
と設定すると新しい時間での温度と濃度 T
十
Cn
+
lが境界条件を満たす.ここで RHS
r
は(
8
.
8
)式の右辺，RHScは(
8
.
9
)式の右辺を差分表示したもの，九， Cbは境界での温度およ
5
6
び濃度である.
境界で，値ではな《法線方向勾配が与えられている場合
l
R
H
.
加
i
)+( 川
T )/&Jd(Xi-Xb)
(
8
.
1
2
)
l
R
H
.
抑
n
J+(C'bn
i-C )/中(
X
i-X
b
)
(
8
.
1
3
)
向 )=
h
(
X
i
)=
n
となる.
以上のように，埋込境界法は運動方程式のみならず，他の輸返方程式にも簡単に適用
できることが分かるが，これまで特に注目して適用された例はみあたらない.次にこの
適用例を示す.
8
.
3浮力噴流の LESシミュレーション
上述のような通常の埋込境界法を一般化した計算法の適用例あるいは検証例として
浮力噴流の計算を行った.これは沿岸域や河川へ放流された下水処理水や発電所冷却水
の拡散などの予測計算の基本形である.ベンチマークとして用いたのは Sharp&Vyas28)
により実験が行われ， Huaie
ta
.
l29)
により 2方程式乱流モデルを用いた
Z
RANS計算が行われている壁噴流で
8
.1に概要が示してあるよ
ある.図 -
f
r
e
es
u
r
f
a
c
e
うにこれは静止流体中の底面に沿
ambient
うように円形噴流が噴射された時
T=EC=O
y
の流れと濃度，温度場の計算である.
放流水は周りの水と異なる密度，温
度を持ち浮力を受けるので，非圧縮
X.
流連続の式と Boussinesq近似の運動
υ
L"
、‘一ーーー一二一一一一参
図8
.1 浮力壁噴流の概要
方程式および，熱拡散，濃度拡散方
程式を空間平均したもの，すなわち
生=一土 VP+丘
二 Pag+V.(vVU+τ
)
+f
Dt
ραρu
¥
7.u=O
与=-v
州
+
v
.(pv円 )+g
5
7
(
8
.
1
4
)
(
8
.
1
5
)
(
8
.
1
6
)
θC
37=-V(uC)+V(
斤 C+J)+h
(
8
.
1
7
)
を解くことになる.ここで u，
，
T Cは空間フィルタ平均した速度，温度，濃度， ρaは周
囲の流体の密度. τ，q，Jはそれぞれサブグリッド応力，サブグリッド熱輸送，及びサ
ブグリッド物質輸送による流束である.これらサブグリッド項には当然モデルが必要で
ある.ここではこれらに標準 Smagorinskyモデル，勾配熱拡散，勾配物質拡散モデルを
用いることにする.詳細は朝倉・中山
3
0
)に記述されているものを用いる.概要は
τ=~ k~~~ d-2v
n~明
3
"
"
5
"
"
"
"
'
6o..J
S
(
8
.
1
8
)
(
8
.
1
9
)
T
V
P一 切
s-phL
一
恥
V
C3
03
E
一
一
T
V
H喧
n
y
-
S=jh+(Vu)T)
(
8
.
2
0
)
Vsgs
J=ysgsVC=-E
了
:
:
.
.
-VC
a
・
(
8.
)
21
sgs
である.ここで dは単位テンソル， Sは速度ベクトル U で定義されるひずみ速度
，
Vs
g
s，
， Ysgsはそれぞれサブグリッド動渦粘性係数，サブグリッド熱拡散係数，サブグリ
ssgs，
ッド物質拡散係数，Pr
s
g
s はサブグリッドプラントル数である
Vs
g
s，k
s
g
sは
Smagorinsky
モデルで算定し， Cs=O.I，Ck=0.09，Pr
0.7と仮定している.
s
g
s=
Vs
g
s= Cs~2ISI
(
8
.
2
2
)
2
}
.
2
1S1
k伊 =Ck/
(
8
.
2
3
)
これらの式の全体的解法は F
r
a
c
t
i
o
n
a
ls
t
e
p法を用い，圧力は HSMAC反復法で解いてい
る.その他の数値計算法の詳細は中山ら
3
1
)
に用いられているものと同じものである.比
較する実験で、は水面はないが，本計算は河川や沿岸での拡散計算に応用するのが目的で
あるので上端は自由水面としている.ただし影響が出ないよう水深は十分大きくとって
いる •
J，g，hの離散化はデルタ関数を Top-hat関数で近似するので，実際は隣接点に分
散することになる.
(
l
)
平坦底面上の壁噴流の計算
まず本手法を，平坦な壁面に接したノズルより周囲の流体と異なる密度 ρjと温度 1
〉
の流体を噴射する時の流れと拡散に適用し，計算法に問題のないことを確認する.噴流
5
8
出口の大きさ，流速，温度等の条件は S
h
a
r
p& Vyasの実験， Huaie
ta
l.の計算に合せ，表
8
.1のような値の計'算ケースを設定する.表で Dt ~はそれぞれ噴流出口の直径，初速
度，九，
ρα は周囲の温度及び密度である • F
dは
ん三円 I~gDj(ρa
-Pj)1Pj
(
8
.
2
4
)
で定義される密度フルード数である. Lは温度希釈率
-立
s
=
二五1
Tj一
九
(
8
.
2
5
)
が 3%となる壁面の噴流出口からの距離 Lで定義される付着距離である.
図8
.2に Case-lの計算結果の瞬時流速と温度分布を示す.計算領域は x方向に予想さ
れる付着距離の約 2倍，横(y)方向に付着距離に相当する領域である.結果は流れと拡
8
.2(
a
)は噴流中心を通る鉛直面内で
散場が静止していない部分のみ図示している.図 の瞬時の流速と温度分布.図 8
.2(
b
)は同じく噴流中心を通る水平面内での分布である.
噴流出口では一様流を設定しているが
流速が遅くなる部分が不安定な塊となり，自然
と乱れが生成される.これらはサーマルの様に上昇している.水平面内の分布図では水
平方向の乱れは噴流全体が浮き上がるまで発達しない.
5
9
40
(
Zm
o
(
a
) 鉛直断面(x=O
)
x(
c
r
r
d
(
b
) 水平断面 (
z
=O
.
5
D
)
図8
.2 流速と噴流濃度分布の瞬時計算結果の例
nHU
ρhu
o 60
T
(
d
e
g
o 20
c
40
60
エ
(
a
) t
=
1
o 60
T
(
d
e
g
.
C
o 20
10
40
X lC
(
b
)t
=
2
図
ー8
・3 流速と噴流濃度分布の
瞬時計算
結果の例(鉛直断面の時間変化)
60
-tE且
nhu
o 60
T
(degハ
o2
0
10
C
40
X
'
60
l
cm
(
c
)t
=
3
o 60
T
(
de
g
.
o20
10
z
40
60
X lC町l
(
d
)t
=
4
図
一8
・3 流速と噴流濃度分布の
瞬時計算結果の例(鉛直断面の時間
変化)
6
2
o 60
T
(
d
e
g.
o 20
10
z
c
40
60
xtcm
(
e
)t
=
5
o6
0
T
(
d
e
g
.
C
)
o20
10
2
C
伽
40
60
x (cm
(
f
)t
=
6
図8
・3 流速と噴流濃
度分布の瞬時計算結果の例(鉛直断
面の時間変化)
6
3
o 60
T
(
d
e
g
.
o 20
G
10
Z
c
40
60
X l
Cπ1
(
g
)t
=7
o60
T
(
d
e
g.
o20
市0
2
C
40
60
Xl
cm
(
h
)t
=8
図
ー8
.
3 流速と噴流濃度分布の瞬時計算結果
の例(鉛直断面の時間変化)
6
4
o 60
T
(
d
e
g
.
C
o 20
10
z
c
60
X lC町E
(
j
)t
=
9
o60
T
(
d
e
g
.
C
o 20
10
c
40
60
X lC町 1
(
j
)t
=
1
0
図8
.3 流速と噴流濃度分布の瞬時
計算結果の例(鉛直断面の時間変化
6
5
)
図8
.4は時間平均した流速と温度分布である.これらの結果を直接比較する実験結果
はないが，出口近傍の壁噴流は渦粘性を仮定した理論値があるので併せてプロットして
いる.この理論は壁極近傍の粘性底層には適用できないので，最大流速の位置より上側
で比較するとよく一致していると言える.温度分布の比較対象はないが，平均時間がや
や短いのかスムーズな等値線になってはいないが，特に問題はない.
壁噴流は Coanda効果で壁に付着する傾向があるが，浮力の効果が大きくなるにつれ
壁から離れていく.この付着距離 Lが計算結果の評価パラメータとなる.これと x=Lで
8
.5に本
の温度希釈率は多くのケースで計測されている.計算も行われているので，図 ta
.
l29)による k
eモデ
計算の Case1，Case2の結果を Sharp& Vyas28)の実験，及び Huaie
ルを用いた計算結果に重ねてプロットしである.本計算結果はこれらのデータと概ね一
致している.
噴流出口の大きさは 4.7mmと流れ場全体に比べると非常に小さいので縦横 4個のセル
で近似しており出口近傍は必ずしも良くない.また壁画近傍の解像度も粘性距離の数倍
と十分ではないが，本研究で用いる数値計算法の基本的に問題ないといえる.
ρ
H
υ
pnu
3
Ver
hof
f
、h
、
、
、
;
一
白、
2
tv
、
、人
1.
5
、、、
﹄伯
王
手
λc
=
i
u
、
、
、
ー
、
、
ー
、、
2.
5
S
.
r
=5D一一一ー
10D
・
LL11
J、
、
¥
tト11 、 1円HI
・-
3
.
5
色
0.
5
O
0
.
2
O
0.
2
0.
4
0
.
6
0.
8
1.
2
U/U
m
(
a
) 出口近傍流速分布
20
30.0
T
Z
の
(deg.
(
cm)
O
20.
0
x (cm)
40
contou:r
・i
ntervαl0.
5 deg.
C
(
b
) 平均温度分布
図8.
4 時間平均流速と温度分布の計算結果，C
a
s
e
l
6
7
240
p
.9
)
ロ EX
‘Huaiω
4
Q ・同
'
一
1
60
80
O
O
20
40
Fa
80
60
(
a
)噴流付着距離
80
60
、、
合
。
、、
q、
*
。
?
。
合
句
40
一
一
一
一
一
一
一
…
'
一
…
.
日
.芝 玄
一
一
…
一
一
一
…
…
.
一
一
一
…
…
一
一
…
.
日
…
…
…
…
.
一
…
…
.
“
…
…
一
…
吋
一
.
一
司
-
一一………一……勾山"一-一…日 百
x
Huai ~ O)
・
・
Case2 z
l
D=70
C制 1 ρ Exp:
)
勾
yJ
xど
壬
支
;
)
(-T恥
h
e
o
α'
1
)
'
.
勺
1
:
'
A
脅 Hua
)
x
'
Tδ
芯
o
σ
'
一
一
…
…
一
一
一
-吾 玄
)
め
9
わ
.
一 T恥
h
e
o
η
'
1
)
'
人
汐
.
川
合
唱
除
凧
20
ロE
Xp.9
ft--p
、事
*
=
z
J
D 35
x
白
一
一
…
一
.
一
一
一
一
一
一
一
一
一
…
…
一
一
…
…
…
一
…
…
一
…
…
.
日
一
玄一
'
Fd
μ
令
(
b
)温度希釈率
8
.5噴流付距離と温度希釈率の実験及び既往計算との比較
図-
6
8
(
2
)
平坦底面に断熱障害物のある場合
次に，埋込境界法の検証のため底面に平たい直方体を置いた場合の計算を行った.こ
=20cm，Lo=6cm，Ho=5cmである.物体は不透過物体である
の物体の位置と大きさは Xo
ので混合物質は通さないが，熱に関しては熱伝導のある場合もありえる.まず熱伝導
i
r
e
c
t
f
o
r
c
i
n
gfを与え，温度場，濃
のない断熱壁の計算を行う.速度の境界には通常の d
度場には熱伝導，物質透過ゼロの Neumann 条件をそれぞれ与える.その他の条件は表
8
.1の Caselと同じである.
8
.6に示す.静止状態から計算を開始し， 6秒後
瞬時速度および温度の計算結果を図 の結果である.物体に衝突することにより噴流が上方に押し上げられる状況が再現され
ている.
6
9
60
T
(
d
e
g.
C
)
20
コ
く
c
:
:
:
>
z
(cm
o
20
x(cm)
40
conlour i
n
t
e
r
v
a
l 27
.
(
a
)鉛直断面
o 60
T
(
d
e
g.
C
)
o20
x(
cm.)
(
b
)水平断面
図8
.6断熱障害物に衝突する浮力噴流の流速と温度分布
7
0
。60
T
(
d
e
g.
o20
o
40
20
(
a
)t
=
1
))
o
40
20
r1
C
r
。
(
b
)t
=
2
図
ー8
.
7断熱障害物に衝突する浮力噴流の流速と温度分布(鉛直断面の時間変化)
7
1
o 60
t~
o
40
20
(
c
)t
=
3
o 60
r
(d
ε
e
o2
0
ト ーてて
o
40
20
(
d
)t
=
4
図8
.7断熱障害物に衝突する浮力噴流の流速と温度分布(鉛直断面の時間変化)
7
2
o6
0
T
(d
e
g
o20
o
、
、
40
20
(
e
)t
=
5
o6
0
T
(d
e
g
o20
L
.
o
40
20
町1
(
f
)t
=
6
図8
.7断熱障害物に衝突する浮力噴流の流速と温度分布(鉛直断面の時間変化)
7
3
o60
T
(
d
e
ご
・
o20
¥
句
、
〈
o
ζ
40
20
(
g
)t=7
o60
T
(
d
e
g.
C
o20
、 『町
、
b
~
o
20
(cm
40
(
h
)t
=
8
図8
.7断熱障害物に衝突する浮力噴流の流速と温度分布(鉛直断面の時間変化)
7
4
ノ
60
T(
deg
o2
0
O
ト
~
:
l
o
40
20
(
i
)t
=
9
b
く
¥
一」三3
o60
T
(
d
e
g・
o20
o
20
40
町1
(
j
)t
=
1
0
図8
.7断熱障害物に衝突する浮力噴流の流速と温度分布(鉛直断面の時間変化)
7
5
、
、
k
(
3
)
平坦底面に定温断障害物がある場合
次に障害物が一定温度を持つ場合，すなわち物体表面で熱伝導があり表面温度を一定
.1の Case1と同じにする.ただし噴
に保つ場合を計算してみる.この場合も条件は表ー8
8
.
1
0
)
流温度 60"Cに対し，物体温度は九=40Cとしている.この場合，温度に対しては式(
0
の九=九に，濃度に対しては式(
8
.
1
3
)で C
'
b
=
O として，境界条件を課することになる.
.8に鉛直断面内および水平断面内での流速と温度分布を示す.流れは
計算結果は図一8
断熱物体の場合と似ているが，温度は物体全面で大きくなりここで浮力が大きくなって
いることが見られる.この場合も詳細を検証する対象がないので定性的判断のみである
が，結果は妥当であるといえる.
もともとプログラミング負荷の小さいところが利点で、あった埋込境界法であるが，効
率等やや犠牲にすることで複雑な状況に対応できる手法であることは他の研究
も示されているが，本提案も有用な一般化になりえる.
7
6
1
6
)3
2
)で
z
(cm
o
20
x(cm)
一
&n
(
a
) 鉛直分布
o
x(
c
r
n
)
(
b
) 水平分布
図8
.8定温障害物に衝突する浮力噴流の流速と温度分布
7
7
o 60
T
(d
e
g.
o20
o
40
(
c
{
a
lt
=
l
o 60
T
(d
e
g
o20
o
40
c
{
b
lt
=
2
図
一8
.9定温障害物に衝突する浮力噴流の流速と温度分布(鉛直分布の時間変化)
7
8
o 60
T
(
d
e
g
.
o20
o
(
cm
40
(
c
)t
=
3
o6
0
T
(
d
e
g.
o20
z
。
40
cm
(
d
)t
=
4
図8
.9定温障害物に衝突する浮力噴流の流速と温度分布(鉛直分布の時間変化)
7
9
o60
T
(
d
e
g
o 20
o
40
(
c
(
e
)t
=
5
o60
T
(
d
e
g
.
C
)
o20
o
40
(
cm
(
f
)t
=
6
図8
.9定温障害物に衝突する浮力噴流の流速と温度分布(鉛直分布の時間変化)
8
0
o6
0
T
(
d
e
g.
o20
z
o
40
(
c
(
g
)t
=
7
o60
T
(
l
c
e
g.
o20
o
内
40
(
c
(
h
)t
=
8
図8
.9定温障害物に衝突する浮力噴流の流速と温度分布(鉛直分布の時間変化)
8
1
o60
T
(
d
e
g.
o20
o
40
(
cm)
(
i
)t
=
9
o6
0
T
(
d
e
g.
o20
z
o
40
(
cm
(
j)
t
=
1
0
図
ー8
.
9定温障害物に衝突する浮力噴流の流速と温度分布(鉛直分布の時間変化)
8
2
第 8章.結論
局所風況予測において，従来のように地形のみをとらえるのではなく，地形上に存在
する構造物・建物の影響を考慮することはとても重要である.本研究では従来から風況
予測の際に用いられてきた境界適合格子法と直交格子を用いて複雑な境界流れに応用
されてきた埋込境界法を融合させ，地形上に存在する構造物を考慮した気流解析手法の
開発と計算手法の検証・応用を行い，さらに埋込境界法を一般化し，単に粘性流の粘着
壁面境界条件のみならず，熱や物質の壁面交換や温度やスカラー値に適用，また直交格
子のみならず，滑らかな境界に適合する座標系に適用し，局所の正確な境界位置で満た
される境界条件と，やや広い領域で境界とみなされる境界を区別する手法を提案してい
る.またその基本的性能や適用法を検討し，気流や熱流放出のシミュレーションへの応
用例も示し以下のような結論を得た.
4章では 2次元丘地形を取り上げ，直交格子を用いて矩形近似法，埋込境界法を用い
て丘地形の境界を表現した.その結果を境界適合格子法の結果と比較することで直交格
子を用いた場合の埋込境界法の計算フログラムの妥当性を示した.
5・6 章においては境界適合格子法に埋込境界法を融合させたプログラムの妥当性を
示すために地形上の任意の位置に存在する任意の大きさの 2次元・ 3次元形状の構造物
境界周りの流れを矩形近似法と比較し
その妥当性を得た.
7章においては気流解析の応用として，実地形上に存在する建物周りの流れの解析を
行い，従来法で、は不可能だ、った境界表現が可能になり，境界適合格子を用いても埋込境
界法の境界表現が有用であるという見識を得た.
実現象をとらえる場合，その流れはレイノルズ数の高い乱流である.その際，境界面
においては粗度を考慮したり，対数則を用いなくては実現象をとらえることは不可能で
ある.しかし，構造物・建物などの境界は流れが直角にぶつかり流れが止まることを考
えると，ノンスリッフ条件を用いることが妥当だと考えられる.
従来の方法では矩形近似法により境界を表現する手法が見られた.矩形近似法は直交
格子を用いることで計算プログラムはとても簡単なもので済むが，地表面においてもノ
ンスリップ境界しか与えることができない.また
高い Re数を対象とする大気流の地
表境界条件に埋込境界法を用いることも同様の理由で適しているとは言えない.
本研究の計算手法は地上に存在する建物の境界はノンスリップ条件を埋込境界法で
表現し，地表面の境界については境界適合曲線座標で表現しているため、地表面の境界
条件に粗度を与えたりすることにより対数則を導入することも可能であり、今後の風況
8
3
予測において実用的なものとなり得る.
8章においては埋込境界法を拡張し，速度がゼ口になる粘着条件，いわゆる D
e
r
i
c
h
l
e
t
条件ではなく，法線方向勾配を指定する Neumann条件への適用法を示している.これ
は熱や物質拡散の地表境界条件に相当し本方法が速度についての境界条件のみでなく
他の輸送式にも使えることを示している.さらに、壁面近傍で噴射される熱噴流に適用
している.この流れは実験データがあり，他の計算例もあり熱を伴う流れの計算検証に
良く用いられる.温度一定の壁や断熱壁の条件が効率よく設定出来ることを示している.
以上のように本論文は，流れの数値解析法に用いる埋込境界法を一般化，拡張し，地
形上気流や熱拡散の乱流シミュレーションに使用可能なものに近づけている.今後気流
のみならず，河川流や温水，淡水の混合する沿岸流のシミュレーション法の改良につな
がると期待される.また本論文で提案する方法を更に一般化すれば，乱流エネルギーや
散逸などの種々の統計量の計算にも応用可能になるものである.
8
4
謝辞
本研究を遂行するにあたり，御指導いただきました神戸大学大学院中山昭彦教授に
は敬意を表すとともに深く御礼申し上げます.また，作業に協力してくれた同講座の
後輩の存在が研究生活の励みになったことは言うまでもありません.著者を支えて下
さった方々に記して感謝の意を表したいと思います.
参考文献
1
) 石原
孟:非線形風況予測モデル MASCOTの開発とその実用化，日本流体力学会，
日本流体力学会誌ながれ， No22，p
p
.
3
8
7
3
9
6，2
0
0
3
.
2
)S
a
n
d
s
t
r
o
r
n
，S
.:WASP-ac
o
r
n
p
a
r
i
s
o
nbetweenr
n
o
d
e
ls
i
r
n
u
l
a
t
i
o
n
sandr
n
e
a
s
u
r
e
r
n
e
n
t
s，Wind
EnergyR
e
p
o
r
tW E9
4
:
2，Dept
.ofMeteorology
，U
p
p
s
a
l
aU
n
i
v
e
r
s
i
t
y
，1
9
9
4
.
.XU，andK
.K
.D
r
o
e
g
e
r
n
e
i
e
r
.:AT
h
e
o
r
e
t
i
c
a
landN
u
r
n
e
r
i
c
a
lStudyofD
e
n
s
i
t
y
3
) Xue，M.，Q
C
u
r
r
e
n
t
si
nN
o
n
c
o
n
s
t
a
n
tS
h
e
a
r Flows，J
o
u
r
n
a
l01TheA
t
m
o
s
p
h
e
r
i
cS
c
i
e
n
c
e
s，Vo.
154，
pp.
1998-2019，1
9
9
7
.
.
，T
a
r
n
u
r
a，T
.
， Okuda，Y，andSanada，S
.:LESoff
l
o
w
so
v
e
ru
r
b
a
n
l
i
k
er
o
u
g
h
n
e
s
s
4
) Nozu，T
b
l
o
c
k
s，The F
o
u
r
t
hI
n
t
e
r
n
a
t
i
o
n
a
lS
y
r
n
p
o
s
i
u
r
n on C
o
r
n
p
u
t
a
t
i
o
n
a
l Wind E
n
g
i
n
e
e
r
i
n
g，
p
p
.
5
7
3
5
7
6，2
0
0
6
.
5
) 村上周三，持田灯，加藤信介，木村敦子:局所風況予測システム LAWEPSの開発
p
.
3
7
5
3
8
6，2
0
0
3
.
と検証，日本流体力学会， 日本流体力学会誌ながれ， No22，p
6
) 内田孝紀，大屋祐二:風況予測シミュレータ RIAM-COMPACT の開発，日本流体
力学会，日本流体力学会誌ながれ， No22，pp.
41
7
4
2
8，2
0
0
3
.
7
) 諏訪好英，藤井修二，湯浅和博，佐野仁美:重合格子法を用いた建築物周辺気流の
数値解析計算アルゴリズムの開発と基本特性の検討，日本建築学会計画系論文集，
488号
， pp.
43，1
9
9
6
.
8
) Guanghao，W.，Binghu，P，andKuroda
，S
.:非圧縮流れの解析に適したカットセル法の
8回数値流体力学シンポジウム， DP-2，2
0
0
4
.
開発，第 1
9
)F
o
r
s
y
t
h
e，G
.E
. andWasow
，W. R
.:F
i
n
i
t
e
D
i
f
f
e
r
e
n
c
eMethods f
o
rP
a
t
t
i
a
lD
i
f
f
e
r
e
n
t
i
a
l
E
q
u
a
t
i
o
n
s，J
o
h
n
-W
i
l
e
ya
n
dSons，I
n
c
.，
p
p
.
1
7
5
2
0
4，1
9
6
0
.
1
0
)中野わかな，畠山望，服部裕司:V
o
l
u
r
n
ePena1
iz
a
t
i
o
n法の誤差低減法，第 25回数値
流体力学シンポジウムラ DI0
・4
，2
0
11
.
I
I
)
P
e
s
k
i
n
，C
.S
.:Flowp
a
t
t
e
r
n
sa
r
o
u
n
dh
e
a
r
tv
a
l
v
e
s，J
comput.Phys，No.10，252， 1
9
7
2
.
1
2
)
P
e
s
k
i
n，C
.S
.:
N
u
r
n
e
r
i
c
a
la
n
a
l
y
s
i
sofbloodf
l
o
wi
nt
h
eh
e
a
r
t，
Jcomput.Phys，No.25，220，
1
9
7
7
.
8
5
，Y
.
，I
s
h
i
h
a
r
a
，S
.，Tanaka
，N
.
:L
a
r
g
e
e
d
d
ys
i
m
u
l
a
t
i
o
no
f
u
n
s
t
e
d
ywindv
e
l
o
c
i
t
yf
i
e
l
d
s
1
3
)H
a
t
t
o
r
i
o
v
e
rah
i
l
lw
i
t
himmersedb
o
u
n
d
a
r
ymethod，TheF
o
u
r
t
hI
n
t
e
r
n
a
t
i
o
n
a
lSymposiumon
p
p
.
3
5
7・360，
2
0
0
6
.
C
o
m
p
u
t
a
t
i
o
n
a
lWindE
n
g
i
n
e
e
r
i
n
g，
1
4
)McQueen，D
.M.，P
e
s
k
i
n，C
.S
.:Shared-memoryp
a
r
a
l
l
e
lv
e
c
t
o
ri
m
p
l
e
m
e
n
t
a
t
i
o
no
ft
h
e
immercedboundarymethodf
o
rt
h
ec
o
m
p
u
t
a
t
i
o
no
fb
l
o
o
df
1
0wi
nt
h
eb
e
a
t
i
n
gmammalian
h
e
a
r
t，JSupercomp，N
o
.
l
l，2
1
3， 1
9
9
7
.
.M.，P
e
s
k
i
n，C
.S
.:
At
h
r
e
e
d
i
m
e
n
t
i
o
n
a
lc
o
m
p
u
t
a
t
i
o
n
a
lmethodf
o
rb
l
o
o
df
10w
1
5
)McQueen，D
C
o
n
t
r
a
c
t
i
l
ef
i
b
e
r
s，J
c
o
m
p
u
t
.
P
h
y
s，No.82，289， 1
9
8
9
.
i
nt
h
eh
e
a
r
t
.I.
，1
. :CombinedI
mmersedB
o
u
n
d
a
r
i
e
s
/
B
S
p
l
i
n
e
sMethodsf
o
rS
i
m
u
l
a
t
i
o
n
so
f
1
6
)Mohd-Yusof
u
a
lR
e
s
e
a
r
c
hB
r
i
e
f
s，
NASAAmes，1
9
9
7
.
Flowi
nComplexG
e
o
m
e
t
r
i
e
s，CTRAnn
. A
.
，
1
7
)F
a
d
l
u
n， E
V
e
r
z
i
c
c
o， R
.， O
r
l
a
n
d
i， P
. and Mohd-Yuso
，f 1
.
Combined
immerced-boundary f
i
n
i
t
e
d
i
f
f
e
r
e
n
c
e methods f
o
rt
h
r
e
e
d
i
m
e
n
s
i
o
n
a
l complex f
10w
.ofComp
t
.P
h
y
s，1
6
1
1，p
p
.
3
5
6
0，2
0
0
0
.
s
i
m
u
l
a
t
i
o
n
s，1
1
8
)下村克弥，尾形陽一:流体ー物体境界上の差分改善による簡明な Immersed-Boundary
法，第 2
5回数値流体力学シンポジウム， DI0-5，2
0
11
.
1
9
)Kim
，J
.，Kim
，D
.a
n
d Choi，H
. :An Immersed-BoundaryF
i
n
i
t
e
V
o
l
u
m
e Method f
o
r
S
i
m
u
l
a
t
i
o
n
sofFlowi
nComplexG
e
o
m
e
t
r
i
e
s，J
o
u
r
n
a
lo
fC
o
m
p
u
t
a
t
i
o
n
a
lP
h
y
s
i
c
s， 1
7
1，
p
p
.
1
3
2・1
5
0，2001
.
2
0
)梶島岳夫:乱流の数値シミュレーション，養賢堂， p
p
.1
l2
1
4
1，1
9
9
9
.
2
1
)
K
a
j
i
s
h
i
m
a，
T
.
， Ohta
，T
.
，O
kazaki，K.a
n
dMiyake，Y
.:H
i
g
h
o
r
d
e
rf
i
n
i
t
e
d
i
f
f
e
r
e
n
c
emethod
10wsu
s
i
n
gc
o
l
l
o
c
a
t
e
dg
r
i
ds
y
s
t
e
m
，T
r
a
n
s
a
ct
.JSMEBVo.
16
3
.No.614，
f
o
ri
n
c
o
m
p
r
e
s
s
i
b
l
ef
p
p
.
3
2
4
7
3
2
5
4
.1
9
9
7
.
.
， Yokojima，S
.andVengadesan
，S
.N
.:C
o
l
o
c
a
t
e
dF
i
n
i
t
eD
i
f
f
e
r
e
n
c
eMethod
22)Nakayama，A
，
i
nG
e
n
e
r
a
lC
u
r
v
i
l
i
n
e
a
rC
o
o
r
d
i
n
a
t
e
sf
o
rS
i
m
u
l
a
t
i
o
nofFlowsOverA
r
b
i
t
r
a
r
yGeometry
C
o
n
s
t
r
u
c
t
i
o
nE
n
g
i
n
e
e
r
i
n
gR
e
s
e
a
r
c
hI
n
s
t
i
t
u
t
eF
o
u
n
d
a
t
i
o
n
，4
3・A，2001
.
2
3
)H
i
r
s
c
h，C.:N
u
m
e
r
i
c
a
lc
o
m
p
u
t
a
t
i
o
nofi
n
t
e
r
n
a
la
n
de
x
t
e
r
n
a
lf
1ows，Volume1
，J
o
h
nW
i
l
e
y
a
n
dS
o
n
s，
NewYork
，1
9
9
8
.
2
4
)Zang，Y
.
，S
t
r
e
e
t
.R
.L
.a
ndK
o
s
e
f
f
，J
.R
.:An
o
n
s
t
a
g
g
e
r
e
dg
r
i
d，
合a
c
t
i
o
n
a
ls
t
e
pmethodf
o
r
t
i
m
e
d
e
p
e
n
d
e
n
ti
n
c
o
m
p
r
e
s
s
i
b
l
eN
a
v
i
e
r
S
t
o
k
e
se
q
u
a
t
i
o
n
si
nc
u
r
v
i
l
i
n
e
a
r
-c
o
o
r
d
i
n
a
t
e
s
.J
C
o
m
p
u
t
.P
h
y
s
.，Vol
.
ll
4，p
p
.
1
8
3
3，1
9
9
4
.
2
5
)宮下康一 :RANSを基にした地形上乱流のシミュレーション，博士論文，神戸大学
大学院自然科学研究科博士課程後期課程建設学専攻， p
p
.
7
1
9，2001
.
2
6
)奥野東，田村哲郎:植生を有する 3次元丘陵地まわりの乱流場の LES解析，日本
風工学誌， Vo
.
13
0，N
o
.
2
(
N
o
.
1
0
3
)，p
p
.
1
4ト1
4
2，2
0
0
5
.
2
7
)川江友勝:シグマ座標を用いた地形上乱流の数値シミュレーション，修士論文，神
8
6
戸大学大学院自然科学研究科博士課程前期課程建設学専攻， 2006
.1
.a
ndVyas，B
.D.: Theb
u
o
y
a
n
tw
a
l
lj
e
t
，P
r
o
c
.I
n
s
t
n
.C
iv
.E
n
g
r
s
.(
L
o
n
d
o
n
)，
2
8
)S
h
a
r
p，1
Vo
.
1
6
3
(
2
)，
p
p
.
5
9
3
6
1
1，1
9
7
7
.
2
9
)Huai，W
x
.，L
i，Z
w
.，Qian，Z
.d
.，Zeng，Y
h
.，Han，，
.
1 Peng，W
q
.:N
u
m
e
r
i
c
a
ls
i
m
u
l
a
t
i
o
n
J
o
u
r
n
a
l01Hydrodynαm
i
c
s
，Vo
.
1
2
2，N
o
.
l，p
p
.
5
8
6
5，2
0
1
0
.
o
fh
o
r
i
z
o
n
t
a
lb
u
o
y
a
n
tw
a
l
lj
e
t，
E
S，土木学会応用力学論文
3
0
)朝 倉 啓 介 ， 中 山 昭 彦 : 乱 流 境 界 層 中 の 混 合 物 拡 散 の L
集
， Vo
.
l1
4，p
p
.
7
5
1
-756，2
0
11
.
3
1)中山昭彦， JeremyD
.B
r
i
c
k
e
r
，本田朔平:沿岸に放流される処理水拡散の観測と LES
シミュレーション，土木学会論文集 B2，Vo
.
1
66，
N
o
.
l，p
p
.
3
5
1
3
5
5，2
0
1
0
.
.a
ndS
o
t
i
r
o
p
o
u
l
o
s，F
. Ah
y
b
r
i
dC
a
r
t
e
s
i
a
n
/
i
m
m
e
r
s
e
dboundarym
e
t
h
o
d
f
o
r
3
2
)G
i
l
n
a
n
o
v，A
s
i
m
u
l
a
t
i
n
gf
l
o
w
sw
i
t
h3D，g
e
o
m
e
t
r
i
c
a
l
l
ycomplex"J C
o
m
p
u
t
.P
h
y
s
. Vo.
1207，457-492，
2
0
0
5
.
8
7