平成１２年度
卒業論文

アイソパラメトリック
リック要素を用いた
流路流れの有限要素解
要素解析

高知工科大学工学部
知能機械システム工学科
知能流体力学研究室

濵田

裕一

目 次

第１章 緒

言

第１節 はじめに１

第２節 本研究の目的と概要

１

第３節 有限要素法の特徴

２

第２章 理論解析
第１節 支配方程式

３第２節 対象領域

４

第３節 境界条件

５

第３章 有限要素方程式の誘導および解法手順
第１節 重み付き残差方程式と弱形式

６

第２節 補間関数

７

第３節 写像関数と形状関数

８

第４節 要素係数マトリックス

１１

第５節 ガウスの求積法

１３

第６節 有限要素方程式

１４

第７節 まとめ

１５

第４章 結果と考察
第１節 要素分割

１６

第２節 解析結果および考察

１７

第５章 結

言

１８

付録 A

１９

参考文献

２７

謝辞

２８

1

第１章

緒 言

第１節 はじめに流体力学の数値解析は，もともと差分法の独壇場であった．一方，有限要素法は固体力学や構造物解析の分野で発展した数値解析法である．有限要素法がはじめて工学上の問題に用いられたのは１９５０年代である．その後，変分法の近似解法へ応用されることによって有限要素法の固体力学・構造物解析への応用は大きな飛躍をとげた．１９６０年代の後半からは，この有限要素法が流体力学の解析手法として認められるようになったが，これは，変分法の拡張としての重み付き残差法を基本とする近似解法が有限要素法と結合されたことが大きく貢献している．特に流体力学における複雑形状物体のまわりの流れ解析には有限要素法が適している．すなわち，有限要素法は汎用性に富み，境界条件の取り込
みが他の手法に比較して容易である．また，有限要素法の基礎は重み付き残差法の弱形式にある．近似解が基本境界条件（ディリクレ境界条件）を満たすこと，および重み関数が基本境界条件上でゼロとなることから，解を容易に求めることができる．

第２節 本研究の目的と概要
以上に述べたように，流体力学にも十分に応用でき，複雑な形状のまわりの流れの解析に適している有限要素法を理解し，特異な性質を有する流体の流れ解析や複雑な形状のまわりの流れ解析をできるようにすることが目的である．２次元の有限要素解析では，解析対象領域の要素分割に三角形もしくは四角形が用いられる．本研究では，アイソパラメト
リック要素を用いた有限要素解析の中で最も基本的な４節点四角形アイソパラメトリッ
ク要素を用いて，２次元におけるポテンシャル流れについて有限要素解析を行う．

2

第３節 有限要素法の特徴
有限要素法の特徴を差分法と比較しながら説明する．両者に共通することは，連続体を細かく分割し多数の独立点における物理量を変数とする連立方程式を解くことにより，離散的に連続体中の物理量を求めようとすることである．差分法では，支配方程式にあらわれる微分量を近似するのに独立点間の物理量の差を用いている．すなわち，微分方程式の未知関数に関するテイラー展開近似を用いて離散化が行われている．一方の有限要素法では，小さな分割された要素によって有限要素方程式を導出している．すなわち，変分法あるいは重み付き残差法を用いる離散化が行われる．また，差分法では原則として直交格子点を採用せざるを得ないのに対し，有限要素法では三角形や四角形などの任意形状を要素として採用することができる．このことにより，有限要素法によれば任意形状の境界形状を表現できることになる．したがって，特に注目したい複雑形状のまわりの流れなども，
その複雑形状部分を細分することによってよりよい近似解を得ることができる．さらに，差分法では取り扱いが不便なノイマン型境界条件（変数の微分値が指定される条件）についても，有限要素法では，自然境界条件をして容易に取り扱うことができる．また，有限要素法では補間関数を任意に選び得ることから，与えられた物理問題に対して最適な選択が可能となっている．
有限要素法は，以上の説明からわかるように一般に汎用性に優れている．しかし，分割する要素数によって計算時間が左右されるという短所もある．一方の差分法は，強い非線形問題（例えば，乱流や衝撃波）には一般に有限要素法よりも計算時間などの点からも有利であると言われている．

3

第２章 理論解析

第１節 支配方程式

一般に２次元の流れにおいて，連続の式は次のように定義される．

¶u ¶v
+
=0
¶x ¶y

(2.1)

ここで u は x 方向速度，v は y 方向速度である．
密度が一定である定常流において，ある関数y(x,y)によって速度が

u=

¶y
¶y
, v=¶y¶x(2.2)

で与えられるとき，連続の式(2.1)は恒等的に満足されていることは容易にわかる．この関数yを流れ関数という．

また，渦度wは次のように定義される．

w=

¶v ¶u
¶x ¶y

(2.3)

ここで，本研究の解析対象は２次元ポテンシャル流れであるので，渦度w=0 である．これを考慮し，流れ関数yを用いて渦度wを表す式を整理することにより，次のような式を得る
ことができる．

¶ 2y ¶ 2y+ 2 =0
¶x 2
¶y

(2.4)

このラプラス方程式(2.4)が本研究の解析対象である２次元ポテンシャル流れにおける支配方程式である．

4

第２節 対象領域
本研究の対象領域は４対１急縮小流路である．寸法は図 2.1 に示し，各頂点に A〜H ま
で記号を設ける．流体は AH から流入し，DE から流出するものとする．
また ABCD，HGFE
は流路壁面であり，一点鎖線は流路中心を示している．

2.0

1.0
B

A

0.25
1.0

y

H

C

D

F

E

G
x

図 2.1

対象領域

図 2.1 に示す対象領域は流路中心に関して上下対称であるため，流れの対称性により実際に計算する解析対象領域は上半分のみとする．したがって寸法は図 2.2 に示すようになり，新たに頂点に J，L，また AB の中点に I，I を通る垂線と中心線の交点に K の記号を
設ける．
2.0
I

A

1.0
B

0.125
0.5

y

C
J

D
L

K
x

図 2.2 解析対象領5

第３節 境界条件
解析対象領域Wの境界は基本境界条件と自然境界条件の２種類に分けられる．基本境界条件は近似関数があらかじめその境界上で満足していることを仮定する境界条件であり，
自由境界条件は近似関数の値を直接与えず，有限要素方程式を誘導する手順において式の中に組み込まれ，近似的に満足される境界条件のことである．
境界条件はポテンシャル分布の絶対的なレベルを定めるために必要な条件で，境界上の少なくとも１点では与える必要がある．ポテンシャル流れの問題は式の偏微分方程式を境
界条件およびのもとに解くことに帰着する．
本解析では，解析対象領域における境界条件はすべて基本境界条件で与えており，図 2.2に示すように流路入口y(y)=y，流路出口y(y)=4y，流路中心y=0，流路壁面y=0.5 である．

入口：y(y)=y
壁面：y=0.5AIB
出口：y(y)=4y
C

J
y

K

流路中心：y=0x

図 2.2 境界条件

D
L

6

第３章 有限要素方程式の誘導および解法手順

第１節 重み付き残差方程式と弱形式
式(2.4)の微分方程式に対して次の積分方程式を考える．
2¶ 2y ù* é¶ yy+òòW êë ¶x 2 ¶y 2 úû dxdy = 0(3.1)

式(3.1)はyが微分方程式の解であれば任意の関数y*=y*(x,y)に対して成立する．言い換えれば，任意のy*に対して積分方程式(3.1)を満足する関数yを求めることができれば，それは微分方程式(2.4)の解である．
式(3.1)をもとにして微分方程式(2.4)の近似解を求めることを考える．関数yが近似解であれば，これは微分方程式を完全には満足しない．
¶ 2y ¶ 2y+ 2 = r ( x, y ) ¹ 0
¶x 2
¶y

(3.2)

この r を残差という．残差が小さいほど，よりよい近似解といえる．式(3.2)より次の式が得られる．

òò y rdxdy = 0
*

W

(3.3)

式(3.3)は，関数y*を重みとして残差が解析領域全体で重み付きで平均的にゼロとなることを意味している．このためy*=y*(x,y)を重み関数と呼ぶ．すなわち，式(3.1)をもとにして近似解を求めるということは，解析領域W内のあらゆる点で残差をゼロにすることはできない代わりに，残差を解析領域W全体で重み付き積分したものがゼロになることを目標として近似解を求める，ということである．このような方法を重み付き残差法といい，式(3.1)
を重み付き残差方程式という．有限要素法は，重み付き残差法に基づいて微分方程式の近似解を求める方法の一つとして位置づけることができ，式(3.1)の近似関数yと重み関数y*に有限要素法独特の関数系を用いる方法である．

7

重み付き残差方程式(3.1)の近似関数と重み関数に対しては一般的に以下の条件が課される．
１．

関数の値が規定されている境界（基本境界条件）上では解の値が分かっているので，近似関数も基本境界上の境界条件を満たしているものとする．

２．

その基本境界上では，重み関数はゼロとなるものとする．

本研究における解析対象領域の場合，領域の両端の境界とも基本境界条件に属するから，その両端で重み関数がゼロであるので，重み付き残差方程式(3.1)を部分積分し整理すると次のように表すことができる．
æ ¶y * ¶y ¶y * ¶y öç÷ dxdy = 0+
W ç ¶x¶x¶y ¶y ÷ø
è

òò

(3.4)

この積分方程式(3.4)を弱形式という．
「弱形式」という用語は，式(3.1)に比べ式(3.4)では関数に対する要求が弱いことに由来する．すなわち，重み付き残差方程式(3.1)の積分が意味あるものとなるには関数が２階微分可能でかつそれが恒等的にゼロにならない必要があるのに対し，式(3.4)では関数が１階微分可能でかつそれが恒等的にゼロにならなければよい．実際，式(3.4)に基づく近似解析のほうが，式(3.1)に基づくよりもより低次の近似関数を用いることができる．

第２節 補間関数
四角形要素を用いて解析対象領域を要素分割した時の１つの要素を代表して，４つの節点(i j k l)からなるものを選び，この要素が占める四角形小領域をWe で表す．４つの節点の座標をそれぞれ<xi yi>,<xj yj>,<xk yk>,<xl yl>とする．また，それぞれの節点におけるyの値をyi,yj,yk,yl で表す．この４つの節点値を補間する関数ye(x,y)は，４つの値を結ぶ曲面を構成する．この曲面の方程式は次のように表すことができる．

y e = a 0 + a1 x + a 2 y + a 3 xy = [1 xyéa 0 ùêa ú
1xy ] ê ú
êa 2 ú
ê ú
ëa3 û

(3.5)

ここで a0,a1,a2,a3 は定数係数である．
この要素補間関数を表す式における定数係数 a0,a1,a2,a3 を求め，式を具体化する．この
手順は次節で写像関係と交えながら説明する．

8

第３節 写像関数と形状関数
有限要素解析に四角形を用いて要素分割した場合，積分方程式を容易に解くために写像という処理を行う．ここで，写像について説明する．解析対象領域を分割する任意形状の四角形要素の補間関数を統一的に記述するために，共通の参照要素を考える．参照要素は辺の長さが２の正方形で，その中心を原点とする(x,h)座標系(以下，正規化座標系)によって要素内の点の位置を記述する．実際の(x,y)座標系(以下，全体座標系)上における四角形
要素内の点(x,y)の座標と正規化座標系上における参照要素内の対応する点(x,h)とを対応づける次の写像関数の具体的な形を求める．

x = x(x , h ), y = y (x , h )x(x,h), y(x,h)

h
1

④

-1

(3.6)

h
③

③

0

④

x

1

x

We
y

①

-1

①

②

②

x
(a)

図3.1

(b)

四角形アイソパラメトリック要素の正規化座標系((a)から全体座標系の要素((b)への写像

図に示したように４つの節点に①①から④④までの番号をふる．これは解析対象領域内のすべての節点を区別するためにつけられる節点番号とは異なり，注目している要素の４つの節点を区別するための局所的な節点番号である．このようにすると実際の要素との節点どうしの対応は明らかである．

9

したがって，式(3.6)の写像関数に対する４つの条件が得られる．
x i = x i (x i , h i ) , y i = y i (x i , h i ) ;(i = 1〜4)(3.7)

ここで(xi,yi)は実際の要素の４つの節点の座標である．また，(xi,hi) は参照要素の４つの節点の座標であり，それぞれ次のように与えられる．
éx 3 ù é 1 ù
éx 4 ù é- 1 ùéx 2 ù é 1 ù
éx 1 ù é- 1 ùêh ú = ê- 1 ú , êh ú = ê- 1 ú , êh ú = ê 1 ú , êh ú = ê 1 úû
ë 4û ë
ë 2û ë
ë 1û ë
ë 3û ë
û
û
û

(3.8)

式(3.7)は x(x,h)と y(x,h)のそれぞれについて４つずつの条件であるから，写像関係式(3.6)
としては４つずつの係数を含んだ関数を考える必要がある．そのような関数として次のものを考える．
x(x , h ) = a 0 + a 1x + a 2h + a 3xh(3.9)

y (x , h ) = b 0 + b 1 x + b 2h + b 3xh
これに式(3.7)，式(3.8)の条件を課すと，８個の係数は次式を満足する必要がある．
é x1 ù é1 - 1 - 1 1 ù éa 0 ù
ê x ú ê1 1 - 1 - 1ú êa ú
ê 2ú = êú ê 1ú ,ê x 3 ú ê1 1 1 1 ú êa 2 ú
ê ú ê
ú ê ú
ë x 4 û ë1 - 1 1 - 1û ëa 3 ûé y1 ù é1 - 1 - 1 1 ù é b 0 ù
ê y ú ê1 1 - 1 - 1ú ê b ú
ê 2ú = êú ê 1ú
ê y 3 ú ê1 1 1 1 ú ê b 2 ú
ê ú ê
ú ê ú
ë y 4 û ë1 - 1 1 - 1û ë b 3 û(3.10)

式(3.10)を解いて式(3.9)に代入すると，写像関係式(3.6)が具体的に次のように表される．
x(x , h ) =44å N (x , h ) x , y(x , h ) = å N (x , h ) y
i

i

i =1

i

i

(3.11)i =1

ここで，Ni(x,h)を形状関数といい，次のように表される．
N 1 (x , h ) =1
1(1 - x )(1 - h ) , N 2 (x , h ) = (1 + x )(1 - h )44(3.12)

1
1
N 3 (x , h ) = (1 + x )(1 + h ) , N 4 (x , h ) = (1 - x )(1 + h )44
また，形状関数 Ni(x,h)は次の関係を満たしている．
N i (x j , h j ) = d ij

(3.13)

ただし，dij はクロネッカのデルタであり次式のように定義される．
é1
ù
ê 1
ú
ú
d ij = ê
ê
1 ú
ê
ú
1û
ë

(3.14)

10

この要素の要素補間関数は，以上のようにして求められた形状関数 NI(x,h)を用いて次のように表される．
4y e (x , h ) = å N i (x , h )y i

(3.15)

i =1

この要素補間関数を全領域にわたって接続したものが結果的に近似関数y(x,y)になる．
このように解析対象領域を要素分割し，近似解を要素ごとに定義された関数によって表現することを離散化という．
式(3.12)に示した要素の形状を表す関数と同じ関数を，写像関係式(3.11)のように要素内の変数の補間にも用いることをアイソパラメトリック写像といい，アイソパラメトリック
写像を用いた有限要素をアイソパラメトリック要素という．要素の形状を表す関数より低次の関数を要素内の変数の補間に用いることはサブパラメトリック写像といい，逆に高次の関数を用いることをスーパーパラメトリック写像という．
補間関数ye(x,h)は辺上では直線であるが要素内では平面ではない．辺上で２節点の値を結ぶ直線分布となることによって，隣接する２つの要素間の補間関数の連続性が保証される．これらの特色は形状関数 Ni(x,h)が式(3.12)に示されているようにxとhの一次式の積になっていることによる．この補間関数を双線形補間関数という．
重み関数に対しても同じ形状関数を用いた関数分布とする．
4y e* (x , h ) = å N i (x , h )y i*(3.16)

i =1

このように近似関数と重み関数とに同じ形の関数を用いる方法をガラーキン法といい，有限要素解析では頻繁に用いられている．

11

第４節 要素係数マトリックス
要素補間関数[式(3.15)]と重み関数[式(3.16)]を弱形式(3.4)に代入することにより，有限要素方程式を導くことができる．要素に関する積分は写像関係を利用して，正規化座標系上で行う．すなわち，被積分関数を一般に G(x,y)とすると，次式のように変数変換して積分
をする．
1

1

òò G( x, y)dxdy = ò ò G( x(x ,h ), y(x ,h )) Jdxdh
We

-1 -1

(3.17)

ここで，J=J(x,h)はヤコビアン
é ¶x
ê ¶x
J (x , h ) = det ê
ê ¶y
êë ¶x
である．式(3.18)の微分¶x ù¶h ú ¶x ¶y ¶x ¶yú=¶y ú ¶x ¶h ¶h ¶x¶h úû(3.18)¶xなどは式(3.11)から具体的に計算することができる．
¶x

ここで弱形式(3.4)を考える．まず，被積分関数を正規化座標系(x,h)の関数に記述することを目標とする．補間関数の微分（¶y e é ¶N 1
=ê
¶x
ë ¶x

¶y e é ¶N 1
=ê
¶y
ë ¶y

¶y e
など）は次のように表される．
¶x

¶N 2
¶x

¶N 2
¶y

¶N 3
¶x

éy 1 ù
ê ú¶N 4 ù êy 2 ú¶x úû êy 3 úê úëy 4 û

¶N 3¶yéy 1 ù
ê ú¶N 4 ù êy 2 úú¶y û êy 3 úê úëy 4 û

ü
ï
ï
ï
ï
ïï
ý
ï
ï
ï
ï
ï
ïþ

(3.19)

さらに，式(3.19)の中の形状関数 Ni の微分は次のように展開される．
¶N i ¶N i ¶x ¶N i ¶h
¶N i ¶N i ¶x ¶N i ¶h
+
=
+
=
,
¶x ¶y ¶h ¶y¶y¶x ¶x ¶h ¶x¶x
ここで，微分¶N i¶Nや i は式(3.12)から容易に求めることができる．
¶x¶h(3.20)

12

一方，

¶x
などは x や y がx,hの非線形関数[式(3.9)]であるため，以下の関係を利用して整
¶x

理する．
é ¶x
ê ¶x
ê
ê ¶h
êë ¶x
重み関数の微分（¶x ù¶y ú 1ú=¶h ú J
¶y úû

é ¶y
ê ¶h
ê
ê- ¶y
êë ¶x

¶x ù
¶h ú
ú
¶x ú
¶x úû

-

(3.21)¶y *など）に対しても同様の展開をすると弱形式(3.4)の被積分関数が¶x

すべてx,hの関数で記述される．
しかしx,hの関数として記述された被積分関数は分母にヤコビアン J がきたりしている
相当に込み入った関数なので，それを積分した関数を求めることは到底不可能である．したがって次節で説明する数値積分法によって係数マトリックスの値を計算する．

13

第５節 ガウスの求積法
アイソパラメトリック要素の係数マトリックスの算出には通常，ガウスの求積法という数値積分法が用いられる．これは次のように関数の積分値を近似するものである．本解析の場合はx方向とh方向の積をとる．
1

1

òò

-1 -1

G (x , h )dxdh @P1

P2

å å G(x ,h )W W
k1 =1 k 2 =1

k1

k2

k1

k2

(3.22)

ここで P は積分点の総数，xk は k 番目の積分点の座標，Wk は k 番目の積分点の重みである．
積分点の個数とその座標および重みの組み合わせを表 3.1 に示す．

積分点数

積分点の座標

重み

0.0

2.0

1

1.0

1

2

±

3

±

3
0.0

8
9

15
5

5
9

表 3.1 積分点の座標と重み

このガウスの求積法によれば積分点の位置における被積分関数の値を算出することによって係数マトリックスを作成することができる．

14

第６節 有限要素方程式
離散化された弱形式(3.4)をこれらにしたがって整理する．ただし，近似関数yと重み関数y*の要素補間関数は
¶N k
¶x
¶N k
¶y

¶N l ù éy i ùú ê ú¶x ú êy j ú¶N l ú êy k ú¶y úû êëy l úû

¶N k
¶x
¶N k
¶y

¶N l ùú¶x ú¶N l ú¶y úûæ ¶y * ¶y ¶y * ¶y öç÷dxdy =+
÷W e ç ¶x
¶
¶
¶
x
y
y
è
ø

é ¶y *
êWe
ë ¶x

é ¶y ù é ¶N iê ¶x ú ê ¶x
ê ¶y ú = êêú ê ¶N iëê ¶y ûú êë ¶y¶N jé ¶y * ù é ¶N iêú êê ¶x* ú = ê ¶x
ê ¶y ú ê ¶N iê ¶y ú ê ¶y
ë
û ë

¶N j

¶x
¶N j
¶y

¶x
¶N j
¶y

ü
ï
ï
ï
ï
ïï
ý
ï
ï
ï
ï
ï
ïþ

éy i* ùê *úêy j ú
êy k* úê *úêëy l úû(3.23)であることを用いる．

òò[= y i* y *j y k* y l*[

= y i* y *j y k* y l*]òòé ¶N iê ¶xêê ¶N jê ¶xêW e ¶N k
ê
ê ¶x
ê ¶N
ê l
ëê ¶x

¶N i ù¶y úú¶N j ú é ¶N iê¶y ú ê ¶xú¶N k ú ê ¶N i¶y ú êë ¶y¶N l úú¶y ûú

é k iiêkê ji
ê k kiêëê k li

k ij
k jj
k kj

k ik
k jk
k kk

k lj

k lk

] òò

é ¶y ùê ¶x úù¶yú ê ¶y ú dxdy¶y û êúëê ¶y ûú
*

¶N j
¶x
¶N j
¶y

¶N k
¶x
¶N k
¶y

k il ù éy i ù
k jl úú êy j úê ú= y*k kl ú êy k úú ê úk ll ûú ëy l û[éy i ù
¶N l ùê úú¶x ú dxdy êy j ú¶N l ú
êy k úê úú¶y ûëy l û

] [k ] [ y ]ee(3.24)e15

ここで，k ij =

1

1 ìæ ¶N¶y¶N i ¶y öæ ¶N j ¶y

¶N j ¶y öò ò íîççè ¶x ¶h - ¶h ¶x ÷÷øççè ¶x ¶h - ¶h ¶x ÷÷ø
i

-1 -1

æ ¶N i ¶x ¶N i ¶x öæ ¶N j ¶x ¶N j ¶x öïü 1
÷ý dxdh÷÷çç ++ çç +¶h ¶x ÷øïþ Jè ¶x ¶h ¶h ¶x øè ¶x ¶h(3.25)

である．
この式(3.25)を用いて要素係数マトリックスを作成したのち，それを全体係数マトリック
スとして組み上げ，境界条件の処理をしてマトリックス方程式を解いて節点速度ポテンシ
ャルの分布を求める．本研究におけるマトリックス方程式の計算は Gauss-Jordan 法で行った．

第７節 まとめ
有限要素法における解析手順を要約すると以下のようになる．
①

連続体の解析対象領域を要素分割し，要素番号，節点番号をつける．本研究では四角形アイソパラメトリック要素を用いる．

②

重み付き残差方程式を作成し弱形式化する．

③

個々の有限要素に着目し，写像関係を利用して有限要素内部の未知変数を要素節点の値を用いて補間する．この近似式には形状関数が用いられる．

④

個々の有限要素に対する離散化した有限要素方程式を作成する．

⑤

要素の係数マトリックスおよび右辺のマトリックスを作成する．本研究における支配方程式では右辺はゼロである．

⑥

全体の係数マトリックスおよび右辺のマトリックスを作成する．これは要素の係数
マトリックスを対応する節点変数成分に重ね合わせることにより行う．

⑦

境界条件の導入と全体の係数マトリックスおよび右辺のマトリックスの変更を行う．

⑧

最終的なマトリックス方程式を解くことにより近似解を得る．

16

第４章 結果と考察

第１節 要素分割
本解析に用いた要素分割を説明する．本研究の解析対象領域は図 2.1 で示したように４対１急縮小流路であり計算領域は上半分である．したがって要素分割も上半分のみとし，
図 4.1 に解析対象領域の要素分割を示す．分割方法は解析対象領域における急縮小部の流路幅を４分割し，その分割幅を１辺とする正方形を用いて解析対象領域を均一的に要素分割している．この時の節点数は 1265，要素数は 1152 である．
A

I

B

C
y
J

D
L

K
x

図 4.1 四角形要素による解析対象領域の要素分割

三角形要素による有限要素解析では，第３章第３節に示している写像関係は用いずに有限要素方程式を誘導するので計算方法に若干違いはあるが，解析対象領域および境界条件が同じである場合，解析結果も本研究とほぼ同等のものが得られるはずである．そこで四角形要素を用いた場合と三角形要素を用いた場合の有限要素解析結果を比較するために，三角形要素を用いた有限要素解析結果を次節に示す．この三角形要素を用いた有限要素解析に用いられた要素分割を図 4.2 に示す．分割方法は図 4.1 に示す各四角形要素に対角線をひいたものであり，この時の節点数は 1265，要素数は 2304 である．
A

I

B

C

D

y
J

K

x

図 4.2 三角形要素による解析対象領域の要素分割

L17

第２節 解析結果および考察
本研究における有限要素解析結果を図 4.3 に示す．また図 4.4 は三角形要素による有限要素解析結果を示す．図 4.3 と図 4.4 はともに急縮小部の流路に着目するため上流部は一部分省略してある．解析対象領域内部に示される線は流線であり，図 4.3 は本研究における４節点四角形アイソパラメトリック要素を用いた有限要素解析，図 4.4 は三角形要素を用いた有限要素解析によって得られた，各節点におけるyの等値線で示されている．また，等値線は流路中心y=0 から値が 0.05 増すごとに表示されている．本研究による有限要素解
析結果は図 4.4 に示す三角形要素による有限要素解析の結果と，ほぼ同等の結果を示している．ともに計算値の詳細は示されていないが，アイソパラメトリック要素を用いた有限要素解析により図 4.3 に示すように２次元ポテンシャル流れにおける妥当な結果を得ることができた．したがって，本研究による解析結果は有効性があるといえる．
I

B

C

D

y
L

K
x

図 4.3 四角形要素を用いた解析結果

I

B

C

D

y
K

L

x

図 4.4 三角形要素を用いた解析結果

18

第５章 結 言本研究では４節点四角形アイソパラメトリック要素を用いた有限要素解析用プログラ
ムを作成し，急縮小流路流れの様子を解析した．また，マトリックス方程式の解法にはGauss-Jordan 法を用いた．三角形要素を用いて解析対象領域を要素分割した場合に比べて，四角形要素を用いて解析対象領域を要素分割した場合は，有限要素方程式の誘導手順に写像関係を有するので若干複雑になる．したがって有限要素解析用プログラムの作成も複雑になるが，図 4.3，図 4.4 に示すように四角形アイソパラメトリック要素を用いた場合と，三角形要素を用いた場合の有限要素解析結果はほぼ同等であるので，本研究において有効性がある有限要素解析プログラムが作成できたとともに，有限要素方程式の誘導，四角形要素を用いて要素分割した場合に特有の写像関係，有限要素解析用プログラムの作成に関する知識などを得ることができた．

19

付録 A

基礎方程式
方程式
流れを２次元非定常流れ，取り扱う流体を非圧縮性粘性流体と仮定したとき，流れの基
礎方程式は次式のようなナビエ・ストークス方程式と非圧縮連続式となる．
¶u i
¶u
¶p
1 ¶ æç ¶u i ¶u j ö÷
+uj i =++¶t¶x j
¶xi Re ¶x j çè ¶x j ¶xi ÷ø

(A-1)

¶u i
=0
¶xi

(A-2)

é x1 ù é x ù
é u1 ù éu ùêu ú = ê v ú , ê x ú = ê y ú
ë 2û ë ûë 2û ë û

(A-3)

ここで，

である．
ここに Re は無次元化のパラメータとしてよく知られているレイノルズ数で，
Re =

rLU
m

(A-4)

である．ただし，rは流体の密度，mは粘性係数，L,U はそれぞれ代表長さと代表流速であ
る．

流速修正法
修正法
流速修正法は，ＳＭＡＣ法に準拠して開発された手法であり，圧力場と速度場を分離して解く分離解法のひとつである．
まず，基礎方程式に対する時間方向の離散化を以下のように行う．
nu in +1 - u in
¶u n
¶p n +1
1 ¶ æç ¶u in ¶u j ö÷
+ u nj i = +
+
Dt
¶x j
¶xi
¶xi ÷ø
Re ¶x j çè ¶x j

¶u in +1
=0
¶xi

(A-5)

(A-6)

20

次に流速の予測子として中間流速 u~i を以下のように定義する．
nìï ¶u n1 ¶ æç ¶u in ¶u j ö÷üïu~i = u in - Dt íu nj i +ý¶xi ÷øïïî ¶x j Re ¶x j çè ¶x j
þ

(A-7)

運動方程式(A-5)の両辺の発散をとり，非圧縮の連続条件式(A-6)を代入して整理すると，次の圧力に関するポアソン方程式が導かれる．
¶ 2 p n +1 1 ¶u~i=
Dt ¶xi¶xi2(A-8)

また式(A-5)から式(A-7)を引くと速度場 u in +1 に関するn +1

¶pu in +1 = u~i - Dt
¶xi

(A-9)

u~i = uˆ i , p n +1 = pˆ(A-10)

が得られる．境界条件はであり，＾は既知量を表す．以上から流速修正法の計算手順は次のようになる．
１．式(A-7)と境界条件式(A-10)より，中間流速 u~i を求める．
２．式(A-8)と境界条件式(A-10)より，圧力 p n +1 を求める．
３．式(A-9)より，速度 u in +1 を求める．

解かなければならない方程式は式(A-7)，式(A-8)，式(A-9)の３つである．これらに重み付き残差法を適用し，まとめると次のようになる．

òW

u i* u~i dW =

n
é
¶u n
¶u i* æç ¶u in ¶u j ö÷ ù1
u i* u in dW - Dt ê u i* u nj i dW+
dW ú
+
W
¶x j
¶xi ÷ø ú
Re W ¶x j çè ¶x j
êW
û
ë

ò

ò

ò

¶u~
¶p * ¶p n +1
1
dW = p * i dW
W ¶x
¶xi
Dt W
¶xi
i

ò
ò

W

u i* u in +1 dW =

ò

ò

W

u i* u~i dW - Dt

ò

W

u i*

¶p n +1dW¶x i(A-11)(A-12)(A-13)

21

補間関数
未知変数を四角形小領域We ごとに表される補間関数を具体化する．四角形小領域内の任意の点の未知関数を４節点での値より補間することを考える．補間関数は次のように表される．

f = N 1f1 + N 2f 2 + N 3f 3 + N 4f 4

(A-14)

ここに，NI は形状関数であり，次のようになる．
N 1 (x , h ) =1
1(1 - x )(1 - h ) , N 2 (x , h ) = (1 + x )(1 - h )44(A-15)

1
1
N 3 (x , h ) = (1 + x )(1 + h ) , N 4 (x , h ) = (1 - x )(1 + h )44
ただし，この形状関数はすでにアイソパラメトリック写像を施しているものである．

有限要素方
要素方程式の誘導
先の重み付き残差方程式(A-11)〜(A-13)において，未知変数であったものは流速 u~i , u i お
よび圧力 p である．そこで四角形有限要素の形状関数を用いて未知変数の補間関数を次のように定める．
u~i = N a u~ai
u i = N uaia(i = 1,2, a = 1〜4)

p = N a pa

ü
ïï
ý
ï
þï

(A-16)

ここで， u~i , u i , p は四角形要素内に分布する未知関数であり， u~ai , uai , pa は４節点における未知関数の値である．次に重み関数 u i* , p * も補間関数を先と同様に定義する．
u i* = N a ua*ip * = N a pa*(i = 1,2, a = 1〜4)

üï
ý
ïþ

(A-17)

補間関数を用いて未知変数と重み関数の近似式が得られたことから，各変数の偏微分の計算があらかじめ行えることになる．

22

例えば，
ü
ï
ï
ý
ï
ï
þ

æ ¶N ¶x ¶N a ¶h öæ ¶N ¶y ¶N a ¶y ö 1¶u ¶N a÷÷ u a÷÷u a = çç a=+u a = çç a¶x¶x
¶h ¶x øè ¶x ¶xè ¶x ¶h ¶h ¶x ø Jæ ¶N ¶x ¶N b ¶h öæ ¶N ¶x ¶N a ¶x ö 1¶v ¶N a÷ va = ç - a=+
va = çç aç ¶x ¶h + ¶h ¶x ÷÷ J va¶h ¶y ÷ø¶y¶yèøè ¶x ¶y(A-18)のようになる (式(3.20)参照）．
以上を用い，四角形要素内で近似された未知変数の補間関数を用いて，先の重み付き残差方程式を整理し，マトリックスの形に書き直すと以下のようになる．
1éù
M ab u~bi = M ab u bni - Dt ê K abgj u bnj +S aibj u bnj úReëûAaibi p bn +1 = -(A-19)

1
H abi u~biDt

(A-20)M ab u bni+1 = M ab u~bi - DtH abi p bn +1(A-21)

ただし，
M ab =

òN a N b dW , K abgj =

We

ò

We

Na N b¶N g¶x j

dWu gni¶N a ¶N b¶N a ¶N b
dW × d ij +dW
W e ¶xW e ¶x¶x ik ¶x kj

S aibj =òAaibi =ò

ò

¶N b
¶N a ¶N b
dW ,H abi =NadW
W e ¶x
We
¶x i
¶x i
i

ò

ü
ï
ï
ï
ï
ý
ï
ï
ï
ï
þ

(A-22)

である．一般に式(A-19)〜(A-21)を有限要素方程式と呼び，式中のマトリックスは係数行列と呼ばれ式(A-22)に示す．

23

この係数行列式(A-22)をさらに整理するとl質量行列 Mab

l

ò

1

1

ò ò N N Jdxdh(A-23)

1 1
æ ¶N g ¶y ¶N g ¶y ö
¶Ng÷dxdhugn
dWugn =
N a N b çç11¶x
¶h ¶x ÷øè ¶x ¶h(A-24)M ab =We

N a N b dW =

-1 -1

a

b

移流項の行列 Kabgj＜i=1,j=1 の場合＞
K abg 1 =

òWe

Na N bòò＜i=1,j=2 の場合＞
K abg 2 =

ò NN
a

We

¶N g

b

¶y

dWu gn =11æ ¶N g ¶x¶N g ¶x öò ò N N ççè - ¶x ¶h + ¶h ¶x ÷÷ødxdhu
a

-1 -1

b

n
g

(A-25)＜i=2,j=1 の場合＞
K abg 1 =

ò NN
a

We

b

1 1
æ ¶N g ¶y ¶N g ¶y ö
¶Ng÷dxdhvgn
dWvgn =
N a N b çç
÷
-1 -1
¶x
¶
¶
¶
¶
x
h
h
x
è
ø

òò

(A-26)＜i=2,j=2 の場合＞
K abg 2 =

ò NN
We

a

b

¶N g
¶y

dWvgn =

1

1

æ ¶N g ¶x¶N g ¶x öò ò N N ççè - ¶x ¶h + ¶h ¶x ÷÷ødxdhv
-1 -1

a

b

n
g

(A-27)

24l

粘性項の行列 Saibj＜i=1,j=1 の場合＞
S a 1b 1 =¶N a ¶N b¶N a ¶N bdW +dW
W e ¶xW e ¶x
¶x
¶x

ò

ò

=2

¶y ¶N a ¶y öæ ¶N b ¶y ¶N b ¶y ö 1
÷ dxdh
çç a
÷ç
-1 -1
¶h ¶x ÷øçè ¶x ¶h¶h ¶x ÷ø Jè ¶x ¶h1

1 æ ¶Nòò(A-28)＜i=1,j=2 の場合＞
1 1 æ ¶N
¶N a ¶N b¶y ¶N a ¶y öæ ¶N a ¶x ¶N a ¶x ö 1
çç a÷ç ÷ dxdh
dW =+W e ¶x-1 -1¶y¶h ¶x ÷øçè ¶x ¶h¶h ¶x ÷ø Jè ¶x ¶hò

S a 1b 2 =òò

(A-29)＜i=2,j=1 の場合＞
Sa 2 b 1 =1 1 æ¶N a ¶N b¶N ¶x ¶N a ¶x öæ ¶N a ¶y ¶N a ¶y ö 1
çç - a÷ç
÷ dxdhdW =+W e ¶y-1 -1¶x¶h ¶x ÷øçè ¶x ¶h¶h ¶x ÷ø Jè ¶x ¶hòòò(A-30)＜i=2,j=2 の場合＞
Sa 2 b 2 =¶N a ¶N b¶N a ¶N bdW +dW
W e ¶yW e ¶y
¶y
¶y

ò

ò

=2

l

1

1 æ

¶N a ¶x¶N a ¶x öæ ¶N b ¶x

¶N b ¶x ö 1

ò ò ççè - ¶x ¶h + ¶h ¶x ÷÷øççè - ¶x ¶h + ¶h ¶x ÷÷ø J dxdh
-1 -1

(A-31)

ポアソン方程式左辺の行列 Aaibi=

1 æ ¶N1¶y¶N a ¶N b¶N a ¶N bdW +dW
W e ¶xW e ¶y¶x¶yòò¶N a ¶y öæ ¶N b ¶y

¶N b ¶y ö 1

Aaibi = Aa 1b 1 + Aa 2 b 2 =ò ò ççè ¶x ¶h - ¶h ¶x ÷÷øççè ¶x ¶h - ¶h ¶x ÷÷ø J dxdh
a

-1 -1

+

1

1 æ

¶N a ¶x¶N a ¶x öæ ¶N b ¶x

¶N b ¶x ö 1

ò ò ççè - ¶x ¶h + ¶h ¶x ÷÷øççè - ¶x ¶h + ¶h ¶x ÷÷ø J dxdh
-1 -1

(A-32)

25lポアソン方程式右辺の行列 Habi
＜i=1 の場合＞
H ab 1 =

ò

Na

ò

Na

We

¶N b
¶x

dW =

æ ¶N b ¶y ¶N b ¶y ö÷dxdhN a çç
-1 -1¶h ¶x ÷øè ¶x ¶hòò(A-33)dW =æ ¶N b ¶x ¶N b ¶x ö÷dxdhN a çç +-1 -1¶h ¶x ÷øè ¶x ¶h(A-34)

11＜i=2 の場合＞
H ab 2 =

We

¶N b
¶y

1

1

òò

となる(式(3.16),(3.19),(3.20),(3.21)参照)．
以上により係数行列の積分が写像後の正規化座標系上で行えるようになる．そこで，先に示した有限要素方程式(A-19)〜(A-21)は最終的に以下のように具体化することができる．
中間流速の式(A-19)＜i=1 の場合＞()(A-35)()(A-36)1ìü
M ab u~b = M ab u bn - Dt í K abg 1u bn + K abg 2 v bn +S a 1b 1u bn + S a 1b 2 v bn ýReîþ

＜i=2 の場合＞
1ìü
M ab v~b = M ab v bn - Dt í K abg 1u bn + K abg 2 v bn +S a 2 b 1u bn + S a 2 b 2 v bn ýReîþ

圧力ポアソン方程式(A-20)(Aa b + Aa b )p b = - D1t (H ab u~b + H ab v~b )
n +1

1 1

2 2

1

2

(A-37)

流速を求める式(A-21)＜i=1 の場合＞
M ab u bn +1 = M ab u~b - DtH ab 1 p bn;1

(A-38)M ab v bn +1 = M ab v~b - DtH ab 2 p bn +1(A-39)＜i=2 の場合＞

26以上に示した式(A-35)〜(A-39)を用い，有限要素解析用プログラムを作成し数値計算を行った．解析対象領域は４対１急縮小流路である(図 2.1 参照)．流速の境界条件として入口の流速 u(y)=-6y2+1.5，v=0，出口の流速 u(y)=-384y2+6，v(y)=0，壁面の流速 u=0，v=0，圧力の境界条件として入口の圧力 p=0 を与えた．マトリックス方程式の計算は本章と同様に
Gauss-Jordan 法を用いた．
この解析の目的はアイソパラメトリック要素を用いた有限要素解析により，計算領域内の速度ベクトル図および圧力分布図を得ることであったが，現段階ではそれに至っていな
い．原因は有限要素解析用プログラムの作成にある．その問題点はいまだ特定できていな
いが，有限要素解析プログラムに含まれている境界条件の入力に関する部分もしくは中間流速を求める式(A-35)および式(A-36)，圧力を求めるポアソン方程式(A-37)，流速を求める
式(A-38)および式(A-39)の入力に関する部分のいずれかに問題があると思われる．

27

参考文献

(1) 日本数値流体力学会有限要素研究委員会，有限要素法による流れのシミュレーション，
シュプリンガー・フェアラーク東京(1998)．
(2) 川原睦人，有限要素法流体解析，日科技連出版社(1985)．
(3) 登坂宣好・大西和榮，偏微分方程式の数値シミュレーション，東京大学出版会(1991)．
(4) 矢川元基・奥田洋司・中林靖，有限要素法流れ解析，朝倉書店(1998)．
(5) 矢川元基，流れと熱伝導の有限要素法入門，培風館(1983)．
(6) 日本機械学会，流れの数値シミュレーション，コロナ社(1988)．
(7) 棚橋隆彦，流れの有限要素解析Ⅰ，朝倉書店(1997)．
(8) 鷲津久一郎ら他３名，有限要素法ハンドブックⅠ基礎編，培風館(1981)．
(9) 鷲津久一郎ら他３名，有限要素法ハンドブックⅡ応用編，培風館(1983)．

28

謝 辞

本研究を遂行するにあたり，終始懇切丁寧にご指導を賜りました蝶野成臣教授，辻知宏
助手に対し，深く感謝いたします．
また高知工科大学知能流体力学研究室の皆さんからは多大なるご援助をいただき，あわせて深く感謝いたします．