1

平成１４年度

修士論文

境界適合法を用いた粘性流動の有限差分解析

高知工科大学大学院
工学研究科

基盤工学専攻知能機械システム工学コース
知能流体力学研究科

国吉

真誓

2

目次

第１章

緒言
１．１

はじめに

1

１．２

格子形成

2１．３

境界適合法

2１．４

研究目的

2

理論および基礎式

3

２．１．１

反復法

3

２．１．２

有限差分法

3

２．１．３

基礎式

4

第２章

理論
２．１

２．２

第３章

6

支配方程式
２．２．１

偏微分方程式

6

２．２．２

座標変換

7

２．２．３

方程式の変換8

２．２．４

計算空間の方程式

10

計算方法
３．１

計算領域

11

３．２

境界条件

13

３．３

計算手順

15

第４章

結果と考察

17

第５章

結言

40

付録

41

参考文献

42

謝辞

43

3

第１章

１．１

緒言

はじめに数値流体力学という学問が始まってから，数値計算における格子形成はそれほど重要視されていなかった．数値流体力学の解法が重要視されるのが当然であり，格子形成は低い評価を与えられていた．理由として，解法に関わる研究は目に見えて結果が出るため研究として認められやすいのに対し，格子形成に関わる研究は結果が表立って見えないため，研究として認められにくかったという事情を指摘することができる．
数値流体力学が普及し始めた頃，研究者・技術者達は，２次元の比較的単純な定番の流れで解析を何度も行う事が多かったため格子形成は容易であった．格子を一度形成すれば，同様のものを何度も利用する事ができたので格子形成がそれほど重要視される事はなく，解法の前に来る付加的な位置づけでしかなかった．
しかしながら，数値流体力学という学問が徐々に発展してくるとともに，その状況も少しずつ変わり，数値流体力学が発展・普及して具体的問題への応用が増加してきた．そして，飛行機や車の周りの空気の流れ，血管内の流れなど複雑な流れ場の高精度なシミュレーションへの要求が高まってきている．流れ場の形状が複雑な場合，離散格子を流れ形状の境界に適合させる必要がある．そのような場合，境界適合法を用いれば複雑な形状の流れに境界が適合し，解析する事が可能になる．また，数値流体力学は，解析事例ごとに異なる形状をして流れ場へどんどん適用されるようになり，格子形成は重要視されるようになった．

１．２

格子形成

流体運動を支配する偏微分方程式を有限差分法や有限要素法で数値計算を行う場
合，解析空間を多数の格子点で離散化する．格子形成はこの解析空間の離散化手法のことを指す．

4

１．３

境界適合法

流れを有限差分法で解析しようとする場合，平行平板管のように境界が直線で表されていれば，直交座標として扱い，簡単に解析ができる．しかしながら，境界が曲線で表されている複雑な流れの場合，解析領域内に格子を作成するために領域を分割すると境界と格子点は一致しない．
その結果，直交座標では解析が困難になり境界付近で計算結果などの精度に問題が出てくる．そこで，境界が曲線で表されている複雑な流れ領域を解析するために，解析する空間を実際の空間とは別に用意し，解析するために境界適合法という手法を用いて解析する．
この境界適合法という手法は，色々な利点が生じる．それは、以下に示す．
・ 現在解析する曲線とは別の曲線で表されている境界でも使用できる．
・ 物体表面が計算空間の境界面に対応しているため，境界条件の与え方が簡単になる．
・ 等間隔格子になるように写像変換すれば差分式および境界条件が簡単になる．

１．４

研究目的

本研究では，解析事例ごとに異なる形状をした流れ場を解析するために現在重要とされている格子形成と数値計算を境界適合法という手法を用いて粘性流動の有限差分解析を行った．

5

第２章

２．１

理論

理論および基礎式

２．１．１

反復法

この手法は，１つの近似値から出発して，方程式を使ってその値を修正しながら正しい解に近づけるという操作を行う．
反復法の長所は，一般に連立１次方程式でなくても適用できること，解に近い値から始めると収束が速いということがあげられる．その反面，方程式の種類や境界条件の種類によっては，収束が非常に遅くなることもあり，あらかじめ何回で収束するか予測がつきにくいなどの欠点がある．

２．１．２

有限差分方程式

有限差分法とは，対象とする領域を格子に分割し，現象を支配している微分方
程式を離散化し，格子点上の値を未知数とする代数連立方程式を解くことで解を求める方法である．
有限差分方程式を以下に示す．

中心差分f ( i +1, j ) − f ( i −1, j )∂f=∂x2∆x(2-1)

f ( i , j +1) − f ( i , j −1)∂f=2∆y∂y(2-2)f ( i +1, j ) − 2 f ( i , j ) + f (i −1, j )
∂2 f
=2
∆x 2
∂x(2-3)f ( i , j +1) − 2 f ( i , j ) + f (i , j −1)∂2 f=2∂y∆y 2(2-4)

6f (i +1, j +1) − f ( i +1, j −1) − f ( i −1, j +1) + f ( i −1, j −1)
∂2 f=
∆x∆y
∂x∂y(2-5)

ただし，式（2-1），（2-2）は１階の差分式，式（2-3），（2-4）は２階の差分式である．

２．１．３

基礎式

２次元の流れでは，連続の式は次のように表される．

∂u ∂v+=0∂x ∂y(2-6)

ここで， u は x 方向速度， v は y 方向速度である．
連続の式が成り立つとき，次の式で定義される関数が存在する．

u=

∂ψ
∂yv=−

∂ψ
∂x(2-7)(2-8)

ここで， ψ は流れ関数と呼ばれる．流れ関数は，恒等的に連続の式を満たす．流れ関数を用いることで，連続の式を考慮に入れる必要がなくなる．
２次元非圧縮性粘性流体の運動方程式は，渦度を用いて以下のように表される．

u∂ω∂ω
1 ∂ 2ω ∂ 2ω+v=+()∂x∂y Re ∂x 2 ∂y 2(2-9)

ただし，ω=∂v ∂u−∂x ∂y(2-10)

ω は，渦度と呼ばれる．
Re =ρULµ(2-11)式（2-7），（2-8）を式（2-9）に代入すると次のように表される．

∂ψ ∂ω ∂ψ ∂ω
1 ∂ 2ω ∂ 2ω+=+()∂y ∂x ∂x ∂y Re ∂x 2 ∂y 2
この式は，渦度輸送方程式と呼ばれる．

(2-12)

7

次に式（2-7），（2-8）を式（2-10）に代入すると次のように表される．

ω=∂∂ψ∂ ∂ψ(−)− ()
∂x∂x
∂y ∂y

∂ 2ψ ∂ 2ψ+= −ω
∂x 2 ∂y 2(2-13)(2-14)

この式は，ポアソン方程式と呼ばれ，流れ関数と渦度との関係を表す．
式（2-12），（2-14）を基礎方程式にとる方法は，流れ関数・渦度法，または，(ψ − ω )法と呼ばれ，２次元非圧縮性粘性流体の数値解を求める標準的な方法の一つとなっている．

２．２

支配方程式

２．２．１

偏微分方程式

境界適合法での座標変換に付随した基礎方程式を誘導しておく．境界適合法を用いた数値計算では、デカルト座標で表されている方程式を曲線座標に変換して差分計算を行う．デカルト座標で表される ( x, y ) 座標を物理空間、曲線座標で表される (ξ ,η ) 座標を計算空間とする．
基礎偏微分方程式は，

f x = ξ x f ξ + η x fη

(2-15)

f y = ξ y f ξ + η y fη

(2-16)

２次元場で考えると写像の変換式は、次のように表すことができる．

ξ x = yη / J(2-17)

ξ x = −x / J(2-18)

η x = − yξ / J(2-19)

η y = xξ / J(2-20)

ただし，J = xξ yη − xη yξ

(2-21)

は変換のヤコビ行列式である．
これらの式（2-17）〜（2-21）を基礎偏微分方程式（2-15），（2-16）に代入する．

8f x = ( yη f ξ − yξ f µ ) / J(2-22)

f y = ( xξ f η − xη f ξ ) / J(2-23)式（2-22），（2-23）により (ξ ,η ) で表される計算空間は，矩形領域となり差分式は単純な式となる．

２．２．２

座標変換

境界適合格子を作成するにあたり，楕円型偏微分方程式を用いた．楕円型偏微分方程式を用いる大きな利点は、解析空間の閉じた境界を指定するとその内部が決まり，物理空間と計算空間の１対１の対応が保証されること，滑らかな格子を作成することができ，境界の不連続性が領域の内部に及ばないことである．物理空間と計算空間の対応を以下に示す．

図１

座標変換

内部格子点の写像関係は楕円型偏微分方程式を満たすことで得られ，ラプラス方程式を用いた．ラプラス方程式は、次のように与えられる．

ξ xx + ξ yy = 0(2-24)

η xx + η yy = 0(2-25)式（2-24），（2-25）を計算空間に変換するために独立変数と従属変数を入れ換

9えると， x, y に関する次の楕円型偏微分方程式を得る．

αxξξ − 2 βxξη + γxηη = 0(2-26)αyξξ − 2 βyξη + γyηη = 0(2-27)

ただし、

α = xη2 + yη2(2-28)

β = xξ xη + yξ yη(2-29)

γ = xξ2 + yξ2(2-30)

である．

２．２．３

方程式の変換１次の変換式を以下に記す．

f x = ( yη f ξ − yξ f µ ) / J(2-31)

f y = ( xξ f η − xη f ξ ) / J(2-32)

ただし，J = xξ yη − xη yξ

(2-33)

は変換のヤコビ行列式である．
２次の変換式を以下に示す．

f xx = { yη ( f x ) ξ − yξ ( f x )η } / J= [ yη {( yη f ξ − yξ f µ ) / J }ξ ] − [ yξ {( yη f ξ − yξ f µ ) / J }η ] / J= ( yη2 f ξξ − 2 yξ yη f ξη + yξ2 f ηη ) / J 2+ {( yη2 yξξ − 2 yξ yη yξη + yξ2 yηη )( xη f ξ − xξ fη )+ ( yη2 xξξ − 2 yξ yη xξη + yξ2 xηη )( yξ f η − yη f ξ )} / J 3(2-34)

f yy = {xξ ( f y )η − xη ( f y )} / J= [ xξ {( xξ fη − xη f ξ ) / J }η ] − [ xη {( xξ f η − xη f ξ ) / J }ξ ] / J= ( xη2 f ξξ − 2 xξ xη f ξη + xξ2 f ηη ) / J 2+ {( xη2 yξξ − 2 xξ xη yξη + xξ2 yηη )( xη f ξ − xξ fη )+ ( xη2 xξξ − 2 xξ xη xξη + xξ2 xηη )( yξ f η − yη f ξ )} / J 3(2-35)

10f xy = {( xξ yη + xη yξ ) f ξη − xξ xη fηη − xη yη f ξξ } / J 2+ {( xξ yηη − xη yξη ) / J 2 + ( xη yη J ξ − xξ yη J η ) / J 3 } f ξ+ {( xη yξξ − xξ yξ ) / J 2 + ( xξ yξη J η − xη yξ J ξ ) / J 3 } f η(2-36)式（2-34），（2-35）を加えると，一般座標でのラプラシアン

f xx + f yy = (αf ξξ − 2 βf ξη + γf ηη ) / J 2+ {(αxξξ − 2 βxξη + γxηη )( yξ fη − yη f ξ )+ {(αyξξ − 2 βyξη + γyηη )( xξ fη − xη f ξ )} / J 3(2-37)

ただし，

α = xη2 + yη2(2-38)

β = xξ xη + yξ yη(2-39)

γ = xξ2 + yξ2(2-40)である．
以上で得られた関係式をデカルト座標で表現された方程式に代入すると計算空間に変換された方程式が求まる．

２．２．４

計算空間の方程式

渦度輸送方程式とポアソン方程式を計算空間での方程式に変換させると次のように表される．

ψ ηω ξ − ψ ξ ωηJ=1[(αω ξξ − 2 βω ξη + γω ηη ) / J 2Re+ {(αxξξ − 2 βxξη + γxηη )( yξ ω η − yη ω ξ )+ {(αyξξ − 2 βyξη + γyηη )( xξ ω η − xη ω ξ )} / J 3 ](2-41)− ω = (αψ ξξ − 2 βψ ξη + γψ ηη ) / J 2 + {(αxξξ − 2 βxξη + γxηη )( yξψ η − yηψ ξ )+ {(αyξξ − 2 βyξη + γyηη )( xξ ψ η − xηψ ξ )} / J 3(2-42)式（2-41），（2-42）を用いて計算を行う．式（2-41），（2-42），境界条件を反復法により解析を行う．

11

第３章

３．１

計算方法

計算領域

本研究で解析を行った領域は，半円状の突起が同 x 位置に存在する場合と異なったx 位置に置かれた場合について解析を行った．また，図の左側が流入，右側が流出とする．流路形状と座標系を図２(a)，(b)にそれぞれ示す．

図２(a)

物理空間

図２(b)

物理空間

12

解析領域は，正方形になっている．計算空間の解析領域の図を以下に示す．

図３３．２

計算空間

境界条件

図３に示すように計算空間は１辺が１の正方形になっており，境界条件は簡単になっている．流入・流出での壁面の境界条件は，変化の大きい場所から十分に距離をとったものと考え、二次元のポアゾイユ流れを与えた．本研究では，物理空間の流入・流出と計算空間の流入･流出の長さを同じにとってあるので速度の方程式を変化させる必要がない．境界条件の式を以下に示す．

流入・流出での境界条件．

u (η ) = −6η 2 + 6η

(3-1)

ψ = ∫ η u (η )dη(3-2)ψ = −2η 3 + 3η

(3-3)

13

ω=−∂u∂ηω = 12η − 6(3-4)(3-5)

上下の壁面に境界条件を有限差分法で与える場合，境界外の点を用いることになる．よって以下の方法で求める．

ψ =1（上部壁面）(3-6)

ψ =0（下部壁面）(3-7)

を与える． ξ 方向速度をψη = 0(3-8)と与え，式（3-6），（3-7），（3-8）と変換した渦度輸送方程式（式（2-12））から上下壁面の境界条件 ω を求める．上下壁面の ω の境界条件を以下に記す．
上部

ω=− 2γψ∆η 2 J 2(3-9)

2γ (1 − ψ )∆η 2 J 2(3-10)

下部

ω=14

境界条件を加えた概略図を以下に記す．

図４(a)物理空間

図４(b)

物理空間

15

次に各領域の格子を示す．

図５図６上下対称に半円がある場合の格子異なる x 座標に半円がある場合の格子

16

３．３

計算手順

１．

計算空間上で物理空間位置を楕円型偏微分方程式により求める．

２．

１．の結果から全格子点上での α ， β ， γ ，および J を計算をする．

３．

図３に示すように計算空間は，１辺が１の正方形となる．この長さは図２(a)，(b)いずれの解析空間で共通であるので１辺１に対する初期値，境界条件を計算し与える．

４．

計算空間上で，式(2-41)および(2-42)を有限差分法を用いて数値計算を行った．

17

第４章

結果と考察

解析領域とレイノルズ数を変化させ，本研究の解析結果を以下に記す．ただし，図は変化部分だけを記す．
まず，各解析領域の流線図を記す．

図 4-1(a)Re=10(b)Re=20(c)Re=25x 方向長さ 10x 方向の半円中心位置 518

図 4-2(a)Re=10(b)Re=20(c)Re=30x 方向長さ 20x 方向の半円中心位置 519

図 4-3(a)Re=10(b)Re=20(c)Re=30x 方向長さ 40x 方向の半円中心位置 1020

図 4-4(a)Re=10(b)Re=20(c)Re=30x 方向長さ 40x 方向の半円中心位置 2021

図 4-5(a)Re=10(b)Re=20(c)Re=30(d)Re=40x 方向長さ 10半円中心位置 5，722

図 4-5(a)Re=10(b)Re=20(c)Re=30(d)Re=40x 方向長さ 20半円中心位置 5，723(a)Re=10(b)Re=20(c)Re=30(d)Re=4024(e)図 4-7Re=50x 方向長さ 40半円中心位置 5，725

次に，各領域の速度ベクトル図を記す．

(a)Re=10(b)Re=2026(c)図 4-7Re=25x 方向長さ 10半円中心位置 527(a)Re=10(b)Re=2028(c)図 4-8Re=30x 方向長さ 20半円中心位置 529(a)Re=10(b)Re=2030(c)図 4-9Re=30x 方向長さ 40半円中心位置 1031(a)Re=10(b)Re=2032(c)
図 4-10Re=30x 方向長さ 40半円中心位置 2033(a)Re=10(b)Re=20(c)Re=3034(d)
図 4-11Re=40x 方向長さ 10半円中心位置 5，735(a)Re=10(b)Re=20(c)Re=3036(d)
図 4-12Re=40x 方向長さ 20半円中心位置 5，737(a)Re=10(b)Re=20(c)Re=3038

図 4-13(d)Re=40(e)Re=50x 方向長さ 40

半円中心位置 5，739

考察・ レイノルズ数が上昇するとともに下流の流れに影響を及ぼす．
・ x 座標に上下対称と異なる x 座標に半円がある場合を比較したとき，上下対称の領域に渦が発生しやすい．
・ 下流の長さが長いほど Re が上昇する．
・ 上下対称の領域の場合，半円の間の領域は図に示すように速度が上昇している．

40

第５章

結言

本研究では，冒頭で述べたように境界適合法を解析に用いた事により境界が曲線で表される複雑な領域の解析を行うことができた．これにより差分法でも複雑な流れ場の解析が可能になる．

41

付録

本研究は，有限差分法で解析が困難な境界が曲線で表される領域の解析が目的である．ここでは，直線格子のみでも解析が可能な直交座標領域を境界適合法で解析したものを付録として示す．境界適合法を用いることで領域が曲線でも直線格子でも解析が可能になる事がわかった．
10.90.80.70.60.50.40.30.20.10345図 A-167格子図

10.90.80.70.60.50.40.30.20.10345図 A-267速度ベクトル図

10.80.60.40.20345図 A-36流線図

742参(1) 河村哲也ら著

考文献

非圧縮性流体解析(1995)東京大学出版(2) 中橋和広・藤井孝蔵著

格子形成方とコンピュータグラフィックス(1995)

東京大学出版(3) 日本機械学会編

流れの数値シミュレーション(4) 高橋亮一・棚町芳弘共著

差分法(1991)(1988)

培風館

日本機械学会43

謝辞

本研究を行うにあたり，終始懇切丁寧なご指導を賜りました蝶野成臣教授，辻知宏
助教授に対し，深く感謝いたします．
知能流体力学研究室の諸氏にも多くの援助をいただき，あわせてここに感謝の意を表します．