C:/大宮司先生関連/大宮司先生原稿作業用1/圧縮性流れの解法３.dvi

(1)

前２章では，まず双曲型微分方程式の特性の理論を説明し，１次元オイラー方程式の初期値問題に関して，特性の理論に基づく３つの流束分離法，流束ベクトル分離法，流束差分分離法，有限体積法について詳述した．次に流束差分分離法と有限体積法を多次元ナヴィエ・ストークス方程式に適用できるように拡張した．本章では更に２つの乱流モデル，

Baldwin-Lomax

の代数モデルまたは

Chien

の２方程式モデルを用いる乱流の計算に拡張する．この拡張によって方程式系の式の数と項の数は増えるが，双曲型の性質は保たれるのでこの拡張は直接的に行うことができる．なお筆者らは，乱流への遷移と乱流境界層の発達はモデルに頼るが，剥離域や後流内の大きい渦は直接計算する手法が必要な精度を有し実用に供し得るものと考えている． 18.1

圧縮性乱流の基礎方程式

圧縮性流れの基礎方程式は質量，運動量，エネルギーの保存の式で次のように表される． @ @t

+

ru u u

= 0

(18.1a)

@ @t u u u

+

r

(

u u uu u u

+

p

) =

r

+

f f f

(18.1b)

@e @t

+

rHuuu

=

r

(

uuu

)

;rqqq

+

fffuuu

(18.1c)

ただしは密度，uuuは流速，e

=

(

+

uuu

2

=

2)

は単位体積当たりの岐点内部エネルギー，pは静圧，H

=

h

+

u u u

2

=

2 = (

e

+

p

)

=は岐点エンタルピー，は粘性応力テンソル，q q qは熱流束，f f fは外力である．式

(18.1)

は圧縮性ナヴィエ・ストークス方程式

(compressible Navier-Stokes equations)

と呼ばれる．

ここで圧縮性流れの数式によく登場する記号について説明する．T は静温度，は比内部エネルギー， h

=

+

p=は比エンタルピーで，完全気体を仮定すれば，p

=

R T d

=

c

v

dT dh

=

c

p

dTである．c

v

は定積比熱，c

p

は定圧比熱，は比熱比，Rは気体定数で，これらの定数の間には

=

c

p

=c

v

R

=

c

p

;c

v

なる関係があり，２つが分かれば他の２つは決まる．c

=

p R Tは音速，M

=

ju u uj=cはマッハ数である．乱流

の直接数値シミュレーション

(DNS, direct numerical simulation)

では，粘性応力テンソル

(viscous stress

tensor)

と熱流束

(heat ux)

qqqの成分は次のようになる．

ij

=

; u

ji

+

u

ij

;

2

3 ij

u

kk

(

ij

= 1

2 3)

(18.2a)

q

i

=

;T

i

=

;

1

;

1

P

r

@c

2

@x

i

(

i

= 1

2 3)

(18.2b)

ただしu

ji

=

@u

j

=@x

i

，

ij

はクロネッカー関数，は粘性係数，は熱伝導率，P

r

=

c

p

=はプラントル数

(Prandtl number)

である．

1

(2)

18.1.1

アンサンブル平均

Navier-Stokes

方程式

乱流の実用解析には，アンサンブル平均ナヴィエ・ストークス方程式

(ensemble-averaged Navier-Stokes

equations)

と乱流モデルの式が用いられる．ナヴィエ・ストークス方程式は乱流にもそのまま適用できる．アンサンブル平均ナヴィエ・ストークス方程式は，ナヴィエ・ストークス方程式の，流れの変数を平均値と乱れ成分に分け，その方程式のアンサンブル平均を取れば導くことができる．ここでアンサンブル平均の意味について考えよう．定常流れの場合には流れ場内のある点の時間平均値と乱れ成分を定義できる．非定常流れの場合にも同様にある点のある瞬間の平均値と乱れ成分を定義したいのであるが，時間平均は取れないので，これに相当のものとしてアンサンブル平均という用語を用いこの場合の平均を定義しているのである．これは頭の中でのことであるが，例えば飛行機が離陸するときに，翼面上の流れのある点の離陸後ある瞬間の流速のアンサンブル平均値は，同じ離陸を何度も繰返し同じ点同じ瞬間の流速を測定し平均を取れば得られるであろう．アンサンブル平均は定常，非定常を問わず用いられが，定常流れの場合には時間平均と同じものになる．圧縮性流れの場合には，密度の重みつきアンサンブル平均，すなわちファーヴル平均

(Favre mean)

が，流速，温度，単位質量当たりで定義される比内部エネルギー，比エンタルピー，エントロピーなどに対して用いられる．例えば密度の一様でない気体の定常乱流を考え，ある点から

1

秒ごとに一定体積の気体をサンプリングしその時間平均温度を出せば，それは単なる時間平均温度ではなく密度の重みつき時間平均温度すなわちファーヴル平均温度になる．上記の量はファーヴル平均を取るのが合理的で，このとき式も簡単になる．ここでファーヴル平均とアンサンブル平均の関係について述べる．乱流の流速，温度等の一つをu で表し，それをアンサンブル平均値と乱れ成分に分ければ u

=

u

+

u 00 u 00

= 0

(18.3)

ただしuはアンサンブル平均，またu 00 は乱れ成分でそのアンサンブル平均はゼロである．またuをファーヴル平均値とその乱れ成分に分ければ u

=

eu

+

u 0 ue

=

u= u 0

=

; 00 u 00 =

(18.4)

ただしeuはファーヴル平均，u 0 はその乱れ成分でそのアンサンブル平均はゼロにはならない

1

．ほかにも u 00

=

00 u 00 u 0

= 0

eu;u

=

00 u 00 = eu

=

ue

=

u

=

u

+

00 u 00 のような関係式が成立する

2

．ナヴィエ・ストークス方程式

(18.1)

のアンサンブル平均を取れば次式が得られる． @ @t

+

ruuu

= 0

(18.5a)

@ @t uuu

+

r

(

uuuuuu

+

p

) =

r

+

fff

(18.5b)

@ @t E

+

r

(

E

+

p

)

u u u

=

r u u u;rq q q

+

f f fu u u

(18.5c)

1 u 0 =u;eu=u; u = 1 (u+u 00 );(+ 00 )(u+u 00 ) = 1 ( u 00 ; 00 u 0 0 )=u 0 0 ; 00 u 00 これよりファーヴル平均の乱れのアンサンブル平均が得られる．なおa =a ab=ab eab= e ab f ab=a e b． 2 eu=(+ 00 ) u =u+ u 0 0 =u

(3)

これらの式の中のアンサンブル平均をファーヴル平均に書き替える． uuu

=

euuu u u uu u u

=

(

e u uue u u u

+2

e u u uu u u 0

+

u u u 0 u u u 0

) =

e u uue u u u

+

g u u u 0 u u u 0 E

=

(

+

u u uu u u=

2) =

e

+12

e u u u e u u u

+12

] u u u 0 u u u 0 e E

(

E

+

p

)

u u u

=

(

h

+

u u uu u u=

2)

u u u

=

e E e u u u

+

p e u u u

+

g h 0 u u u 0

+

g u u u 0 u u u 0 e u u u

+12

^ u u u 0 u u u 0 u u u 0 これらの関係を用いれば，次のアンサンブル平均ナヴィエ・ストークス方程式

(ensemble-averaged

Navier-Stokes equations)

が導出できる． @ @t

+

r e u u u

= 0

(18.6a)

@ @t uuue

+

r

(

uueueuuu

+

p

) =

r

(

f ; g uuu 0 uuu 0

) +

f ff

(18.6b)

@ @t e E

+

r

(

e E

+

p

)

uueu

=

r ;

(

f ; g u uu 0 u uu 0

)

euuu;qqq; g h 0 u u u 0

+

fffuuu

(18.6c)

なお上式の導出に際しては３重相関項 ^ u u u 0 uuu 0 uuu 0 =

2

は省略された

3 .

ところで乱流の計算では，アンサンブル平均ナヴィエ・ストークス方程式

(18.5)

ないし

(18.6)

を用いず通常のナヴィエ・ストークス方程式

(18.1)

を用い，その拡散項すなわち応力テンソルの成分

ij

と熱流束の成分q

i

を

ij

=

; u

ji

+

u

ij

;

2

3 ij

u

kk

;u 0

i

u 0

j

(

i j

= 1

2 3)

q

i

=

;

1

;

1

P

r

c

2 _i

₊

h 0 u 0

i

(

i

= 1

2 3)

のように置くことが多い．ただし;u 0

i

u 0

j

=

; g u 0

i

u 0

j

はレイノルズ応力

(Reynolds stress)

，h 0 u 0

i

=

g h 0 u 0

i

は乱流熱流束

(turbulent heat ux)

の成分で，共に乱流の乱れによる拡散項である．薄剪断層の層流域では

右辺第

1

項の分子拡散のみが働き，乱流域では第

2

項の乱流拡散が支配的である．次にアンサンブル平均ナヴィエ・ストークス方程式ではなく通常のナヴィエ・ストークス方程式をそのまま用いることの意味について説明する．もとのナヴィエ・ストークス方程式

(18.1)

は，

ij

とq

i

を上式のように置き，更に ) uuu)uuue e) e E p)p H ) e E

+

p= と置けば，ファーヴル平均を用いたアンサンブル平均ナヴィエ・ストークス方程式

(18.6)

になる．乱流モデルを用いる圧縮性乱流の計算では，表面的には通常のナヴィエ・ストークス方程式を解いているように見えるが，実はアンサンブル平均ナヴィエ・ストークス方程式を解いているということである．したがって乱流の計算における式

(18.1)

の例えばeは，単位体積当たりの岐点内部エネルギーではなく e) e E

=

; e

+12

uueueuuu

+12

uuu 0 uuu 0 3 乱流では平均流に対して乱れ成分は必ずしも小さくなくこの項は無視できないのであるが，そのモデル化はことを煩雑にするので，拡散を表す項はレイノルズ応力項に代表させこの項を省略したということである．なおこれに関しては筆者の誤解があるやも. e h= ; e + f p= f p== ; p= =p

(4)

ということになる．上記のように記号を対応づけた場合に，アンサンブル平均した状態方程式は p

=

R e T de

=

c

v

d e T d e h

=

c

p

d e T

(18.7a)

また熱力学的状態量の関係式等は p

= ~

e

= ~

e E; f

2 (18.7b)

e H

=

e E

+

p

=

e E; f

2 = 1~

(

f c

2 +

f

2 )

(18.7c)

のようになり，形式的に通常のものと同じになる．ただし

~

=

;

1

f

2 = 12~

(

uueuuueu

+

] u uu 0 uuu 0

)

f c

2 =

R e T ec

=

q f c

2

これから先は平均を表す

` '

と

`

e

'

は省略することにする．

18.1.2

座標成分で書かれた圧縮性

Navier-Stokes

方程式

圧縮性ナヴィエ・ストークス方程式

(18.1)

はデカルト座標系

(cartesian coordinates)

xxxでは次のようになる． @q @t

+

@F

i

@x

i

+

D

+

g

= 0

(18.8a)

q

=

0 B B B B B B B @ u

1

u

2

u

3

e 1 C C C C C C C A F

i

=

0 B B B B B B B @ u

i

u

1

u

i

+

1 i

p u

2

u

i

+

2 i

p u

3

u

i

+

3 i

p Hu

i

1 C C C C C C C A

(

i

= 1

2 3)

(18.8b)

D

=

; @ @x

i

0 B B B B B B B @

0

1 i

2 i

3 i

ij

u

j

;q

i

1 C C C C C C C A g

=

; 0 B B B B B B B @

0

f

1

f

2

f

3

f

i

u

i

1 C C C C C C C A ただし上式のqは未知変数，F

i

は流束，Dは拡散項，gは外力項である．

次に一般曲線座標系

(general curvilinear coordinates)

の場合を考える．式

(18.1)

にヤコビアンJを乗じ，Jr

=

@ @

i

J

(

r

i

)

なる関係を用い，反変速度U

i

=

uuur

i

を導入すれば @ @t J 0 B B @ uuu e 1 C C A

+

@ @

i

J 0 B B @ U

i

uuuU

i

+

r

i

p HU

i

1 C C A

=

@ @

i

(

Jr

i

)

0 B B @

0

u u u; q q q 1 C C A

+

J 0 B B @

0

g gg g g gu u u 1 C C A

(5)

これより次の一般曲線座標系の圧縮性ナヴィエ・ストークス方程式が導かれる． @q

^

@t

+

@F

^

i

@

i

+ ^

D

+^

g

= 0

(18.9a)

^

q

=

J 0 B B B B B B B @ u

1

u

2

u

3 ^

i

=

J 0 B B B B B B B @ U

i

u

1

U

i

+

i

1

p u

2

U

i

+

i

2

p u

3

U

i

+

i

3

p HU

i

1 C C C C C C C A

(

i

= 1

2 3)

(18.9b)

^

D

=

; @ @

i

J

ij

0 B B B B B B B @

0 j

1 j

2 j

3 jk

u

k

;q

j

1 C C C C C C C A g

^

=

;J 0 B B B B B B B @

0

f

1

f

2

f

3

f

k

u

k

1 C C C C C C C A ただしq

^

は未知変数，F

^

i

は流束，D

^

は拡散項，

^

gは外力項である

4

．またJ

=

@

(

xyz

)

=@

(

)

，

ij

=

@

i

=@x

j

は変換のヤコビアンと測度である．式

(18.9)

は，式

(18.8)

のx x xの微分をの微分に変え，曲線座標格子を用いる数値計算に適するようにしたものである．ここで式

(18.9)

の意味について少々説明する．数値解析では，物理空間内の曲線座標格子

(curvilinear

coordinate grid)

は通常格子間隔

1

の立方体格子に写像され，写像空間内の格子セルに対して差分法または有限体積法による計算が行われる．まずヤコビアンJの意味を考えよう．図

18.1

を参照すれば， J

=

x

y

z

(

xxx

100

; xxx

000 )

(

xxx

010

; xxx

000 )

(

xxx

001

; xxx

000 )

したがってJは空間内の単位立方体要素のx x x空間における体積，すなわち物理空間の局所体積：写像空間の相当体積の比を表している．なおここではx

1 =

x x

2 =

y x

3 =

z

1 =

2 =

3 =

としている．次にpにかかる係数，例えばJr の意味を考えよう．下図を参照すれば， Jr

=

J

x

J

y

J

z

=

y

z

x

y

!

(

x x x

010

; x x x

000 )

(

x x x

001

; x x x

000 )

H H j : 6 H H H H H HH H H H H H H H H H H H H H 000 100 010 001 H H j U000 H H j U100 : V 000 : V010 6 W000 6 W 001 図

18.1:

空間内の単位立方体格子 4 デカルト座標系から一般曲線座標系への変換に関して不明の点は17.2.2項を参照ください．

(6)

Jr は，空間内の単位立方体要素の

=const.

面に相当するxxx空間内の面の面積ベクトルである．つまりその大きさは物理空間の

=const.

面上の局所面積：写像空間の相当面積の比を表している．式

(18.9)

は空間の単位立方体要素の流れに対して解釈すれば次のようになる．その密度，物理空間の運動量，単位体積当たりの岐点内部エネルギーは物理空間のものにヤコビアンJを乗じたものになり，それらの流束は反変速度U U Uで要素へ流入または流出する．なお反変速度は空間内の流れの速度である．また要素境界

i

=const.

に作用する圧力，粘性応力，この面を通る熱流束は物理空間のものにJjr

i

jを乗じたものになる．xxx空間と空間の対応する２つの要素に関して，要素内の流れの質量，岐点内部エネルギー，物体力は同じになる．また対応する境界面に作用する力の成分と通過する熱量も同じになる．次に空間の流れの運動方程式を作る．反変速度成分U

l

は流速uuuにr

l

を乗じたものであるから，この運動方程式は式

(18.1b)

にJr

l

を乗じることによって導かれる．なおその導出過程では Jr

l

uuu

=

JU

l

Jr

l

(

ru u uu u u

+

p

) =

r

l

@ @

i

J

(

U

i

u u u

+

r

i

p

) =

@ @

i

J

(

U

i

U

l

+

g

il

p

)

;J

(

U

i

u u u

+

r

i

p

)

@ @

i

r

l

のような計算が行われる．上記の質量とエネルギーの式にこの運動方程式を加えた写像空間の流れのナヴィエ・ストークス方程式は次のようになる． B @q

^

@t

+

@F

^

i

@

i

+ ^

D

+^

g

=

@q

~

@t

+

@F

~

i

@

i

+ ~

R

+ ~

D

+~

g

= 0

(18.10a)

~

q

=

1

U

2

U

3 ~

i

=

i

U

1

U

i

+

g

i

1

p U

2

U

i

+

g

i

2

p U

3

U

i

+

g

i

3

p HU

i

1 C C C C C C C A

(

i

= 1

2 3)

~

R

=

;J

(

U

i

uuu

+

r

i

p

)

@ @

i

0 B B B B B B B @

0

r

1

r

2

r

3

0

1 C C C C C C C A

(18.10b)

~

D

=

;J @

j

@x

k

@ @

j

ik

0 B B B B B B B @

0

1 i

0

2 i

0

3 i

0

u

i

1 C C C C C C C A ; @ @

i

Jg

ij

k @T @

j

0 B B B B B B B @

0

1

1 C C C C C C C A g

~

=

;Jf f f 0 B B B B B B B @

0

r

1

r

2

r

3

u uu 1 C C C C C C C A B

=

0 B B B B B B B @

1 0

0 0 0

0

1

2

1

3

0

2

1

2

3

0

3

1

3

2

3

0 0 0

0 0 1

1 C C C C C C C A

(18.10c)

(7)

ただし式

(18.10b)

のq

~

は未知変数，F

~

i

は流束，R

~

は付加項，D

~

は拡散項，g

~

は外力項で，式

(18.10c)

のB は流速u u uを反変速度U U Uに変換する行列である．またg

ij

=

ik

jk

は反変測度テンソルの成分で，添字

0

の付いている量は微分しないものとする．@

(

a

0

b

)

=@

=

a

0

@b=@ ．

式

(18.10)

の運動方程式は質量流束

(mass ux)

JUUUの運動方程式である．この質量流束の

i

成分JU

i

は

空間の単位立方体格子をx x x空間に写像した格子セルの

i

= const.

面を通る質量流量である．この運動方程式の

3

成分の式は，式

(18.9)

の保存形運動方程式の

3

成分の線形結合を取ったものであるが非保存形である．式

(18.10)

は，対流項を保存形に書換えたもので，そのために付加項R

~

が生じている．この付加項は無視できないがその値は滑らかな格子では小さい．

18.1.3

乱流モデル

|Baldwin-Lomax

の代数モデルと

Chien

の

2

方程式モデル

乱流モデルの選び方とその適用限界については第10章を参照されたい．

Baldwin-Lomax

モデルは，代数モデルのひとつで，この種のモデルでは乱流を特徴づける変数，このモデルでは混合長lが，その断面内の状態

(

流速プロフィル，圧力勾配，壁面曲率など

)

だけから決まる．これに対し

Chien

モデルは低レイノルズ数型k;"モデルのひとつで，この種のモデルでは乱流を特徴づける乱れの変数，運動エネルギーkと散逸率"の輸送方程式が解かれ，代数モデルとは異なり上流側の影響が直接計算結果に反映される．しかしながら渦粘性近似を置いているため，レイノルズ応力の１つの成分が卓越している単純な流れには良いが，複雑乱流には不十分である．このような流れには代数応力モデルなどのより高度のモデルを選ばなければならない．低レイノルズ数型モデルは壁面近傍の低レイノルズ数域にも適用でき，壁法則を補う必要のないモデルである．ここには，

Baldwin-Lomax(1978)

の代数モデルと

Chien(1982)

の

2

方程式モデルを示すことにする．これらの渦粘性型モデルでは応力テンソルと熱流束q q qの成分は次のように置かれる．

ij

= (

+

t

)

u

ji

+

u

ij

;

2

3 ij

u

kk

;

2

3 ij

k

(

ij

= 1

2 3)

(18.11a)

q

i

=

;

1

;

1

P

r

+

t

P

rt

@c

2

@x

i

(

i

= 1

2 3)

(18.11b)

ただし，

t

は渦粘性係数

(eddy viscosity)

，k

=

u 0

i

u 0

i

=

2

は単位質量当たりの乱れの運動エネルギー

(turbulent

kinetic energy)

，また P

rt

は乱流プラントル数

(turbulent Prandtl number)

で，空気に対してはP

r

=

0

:

72

P

rt

= 0

:

9(

境界層

)

= 0

:

5(

自由剪断流

)

である．

式

(18.11a)

と

(18.11b)

は，それぞれ層流の場合の式

(18.2a)

または

(18.2b)

に同形の乱流の項を加えた

もので，層流のところでは

t

= 0

で乱流項の寄与はゼロになる．式

(18.11a)

のレイノルズ応力

(Reynolds

stress)

項は，レイノルズ応力の非等方テンソル成分が平均歪み速度の非等方テンソル成分に比例するものと

してそれらの成分に渦粘性係数

t

を乗じたもので，速度勾配拡散の式である．また式

(18.11b)

の乱流熱流

束

(turbulent heat ux)

項は，平均温度勾配;@T=@x

i

に渦熱伝導率

(thermal eddy viscosity)

t

=

c

p

t

=P

rt

を乗じたもので温度勾配拡散の式である．これらの勾配拡散モデルでは，すべてが拡散係数

t

のモデル化

にかかっている．そのモデル化は発達段階の乱流境界層や

2

次元ジェット流のようにレイノルズ応力の１つ

の成分が卓越している場合には各種の効果を含めて可能であるが，これらの薄剪断層内に横流や大きい縦渦のある場合には難しくなる．

(8)

Baldwin-Lomax

モデルは

2

層代数渦粘性モデルで，渦粘性係数は次のように置かれる．

t

=

(

t

)

inner

(

dd

crossover

)

(

t

)

outer

(

d>d

crossover

)

(18.12a)

ただしdは壁面からの距離，d

crossover

は

(

t

)

inner

= (

t

)

outer

になるところのdである．内層には

Prandtl-van Driest

の式が用いられる．

(

t

)

inner

=

l

2

jj

(18.12b)

ただしl

=

df

1

;

exp(

;y

+

=A

+

)

gは混合長，

=

ru u uは渦度である．

= 0

:

4

は

Karman

定数，A

+

= 26

， y

+

=

u

d=

w

=

p

w

d=

w

は壁面からの無次元距離，u

=

p

w

=

w

は摩擦速度，

w

は壁面剪断応力である．他方外層には

Clauser

の式の代わりに次式が用いられる．

(

t

)

outer

= 1

:

6

KF

wake

F

Kleb

(

d

)

(18.12c)

ただしK

= 0

:

0168

は

Clauser

定数である．また

F

wake

= min

;

d

max

F

max

0

:

25

d

max

u

dif2

=F

max

F

Kleb

(

d

) =

h

1+5

:

5

0

:

3

d d

max

6

i ;

1

d

max

はF

(

d

) =

djjf

1

;

exp(

;y

+

=A

+

)

gが薄剪断層を横断して最大値F

max

になるところのdである．ただ

し後流では

exp(

;y

+

=A

+

) = 0

とする．またu

dif

=

juuuj

max

;juuuj

min

は薄剪断層を横断しての流速の大きさの

差で，境界層ではju u

uj

min

= 0

とする．F

Kleb

は

Klebano

の間歇係数である．遷移については，薄剪断層を

横断して

(

t

)

max

<

14

1

(18.13)

ならば，その位置では薄剪断層内の流れは層流で薄剪断層を横断して

t

= 0

とする．

Baldwin-Lomax

モデルを用いる場合には，前章に述べた解法をほとんどそのまま使用できる．ただし境界に沿って

(

t

)

max

を計算し条件

(18.13)

から遷移点を求めることが必要になる．またその下流域では境界付近の流れが乱流になるので，ナヴィエ・ストークス方程式の拡散項は渦粘性近似の式

(18.11)

と

Baldwin-Lomax

モデルの式

(18.12)

を用い計算される．

Chien

モデルは低レイノルズ数型k;"モデルで，このモデルはもともとy

= const.

壁面に沿う非圧縮性流れに対して開発されたものであるが，

Chien

モデルでは一般曲線座標系の圧縮性流れに適用できるよう次のように書換えられている．

t

=

c

f

k

2

="

(18.14a)

乱れの運動エネルギーkとその散逸率"は次の輸送方程式を解くことによって求められる． @q

^

t

@t

+

@F

^

ti

@

i

+ ^

D

t

+^

g

t

= 0

(18.14b)

(9)

ただし

^

q

t

=

J k " ! F

^

ti

=

JU

i

k " ! D

^

t

=

; @ @

i

Jg

ij

0 B B @

(

+

t

)

@k @

j

+

t

@" @

j

1 C C A

(18.14c)

^

g

t

=

;J 0 @ P;";D " k

(

c

1

P;c

2

f";f

3

D

)

1 A P

=

;u 0

i

u 0

j

@u

i

=@x

j

は乱れの発生項，;u 0

i

u 0

j

はレイノルズ応力，D

= 2

k=d

2

は壁面におけるエネルギー散逸率，dは壁面からの距離である．またモデル定数とモデル関数は次のように与えられる． c

= 0

:

09

c

1 = 1

:

35

c

2 = 1

:

8 = 1

:

3

f

= 1

;e ;

c

3

y

+ f

= 1

;

0

:

4

1

:

8

e ;

(

Re

t

=

6)

2 f

3 =

e ;

c

4

y

+ c

3 = 0

:

0115

c

4 = 0

:

5

， R e

t

=

k

2

=" R e

t

は乱流レイノルズ数である．このモデルは壁近傍の低レイノルズ数流れにまで適用できるが，壁面の隣接格子点は粘性底層内のy

+

du

=

= 1

4

に取らなければならず，壁面近くに細かい格子が要求されることになる．実際問題として壁面近くの格子は特別の注意を払わなければ非常に細長いものになり，このことが壁面に直立する衝撃波の捕獲に悪影響を及ぼすようである．また通常上流境界にkと"の小さい初期値が適宜与えられる．この主流の乱れは境界層内の流れに外乱となって作用しその流れを不安定にしバイパス遷移を惹き起す．なお上流境界にkと"の初期値を与えることは必要なことで，プログラムにこのモデルを組込んだだけでは流れの全域がk

=

"

= 0

になり乱流の計算にはならない． k;"モデルは，良く検証され信頼性の高いモデルであるが，kは乱れの特性速度qの

2

乗，"はqの

3

乗に比例する量で，流れが大きく乱れるところではこれらの量は極端に大きくなり，計算が不安定になるという難点を抱えている．安定化のためには，層流として計算を始める k方程式の発生項Pの大きさを散逸項と対比しながら抑える kと"の輸送方程式をナヴィエ・ストークス方程式と同時に解くなどの対策を講じることが必要である．なお計算格子を細かくすることはkと"の最大値を更に大きくし問題解決にはならない．上記の第

3

の対策のためには，ナヴィエ・ストークス方程式にkと"の輸送方程式が加わった場合の流束差分分離法，有限体積法，近似因子法の諸式を導出することが必要となる．

(10)

18.2

線形化と対角化

以下には低レイノルズ数型k;"モデルを用いるアンサンブル平均圧縮性ナヴィエ・ストークス方程式の解法について述べる．レイノルズ応力輸送方程式モデルなど，他の輸送方程式モデルを用いる場合については省略するが，これらの場合も同様の解法で解を求めることができる．以下本章にはkと"の方程式が加わった場合に，前章に述べた圧縮性ナヴィエ・ストークス方程式の２つの解法，流束差分分離法と有限体積法，それから計算効率を上げる方法がどのようになるのかについて一通り説明する．本節には２つの解法の説明に入る前の下準備として，

1

次元オイラー方程式にkと"の方程式が加わった場合の線形化

(linearization)

，非保存形化，対角化

(diagonalization)

について説明する．

18.2.1

デカルト座標系の方程式

前章で線形化・対角化の対象として用いられたオイラー方程式はナヴィエ・ストークス方程式から拡散項と外力項を除外したものであったが，これに加えられるべき相当のkと"の方程式は前出のkと"の輸送方程式

(18.14)

から拡散，発生，散逸の項を除外した輸送項のみの方程式である．これらの方程式は纏めて書けば次のようになる． @ @t

+

ruuu

= 0

(18.15a)

@ @t u u u

+

r

(

u u uu u u

+

p

) = 0

(18.15b)

@e @t

+

rHuuu

= 0

(18.15c)

@k @t

+

rkuuu

= 0

(18.15d)

@" @t

+

r"u u u

= 0

(18.15e)

なおこれらの方程式系は双曲型の性質を保持している．この方程式系から１つの空間次元の部分のみを取り出せば次式が得られる． @q @t

+

@F

i

@x

i

= 0

(

i

= 1

2 3)

(18.16a)

ただし添字iは

1

または

2

または

3

で，ここではアインシュタイン総和規約には従わないものとする．未知変数ベクトルqと流束F

i

は次のようになる． q

=

0 B B B B B B B B B B B B B @ u

1

u

2

u

3

e k " 1 C C C C C C C C C C C C C A F

i

=

i

u

1

u

i

+

1 i

p u

2

u

i

+

2 i

p u

3

u

i

+

3 i

p Hu

i

ku

i

"u

i

1 C C C C C C C C C C C C C A

(

i

= 1

2 3)

(18.16b)

最後の

2

行が乱流モデルのものである．

(11)

流束F

i

はqの成分の

1

次の同次式であるから，オイラーの同次関係を用いれば dF

i

=

@F

i

@q dq

=

A

i

dq

(18.17)

F

i

=

A

i

q dA

i

q

= 0

となり，式

(18.16)

は次のように線形化することができる． @q @t

+

A

i

@q @x

i

= 0

(

i

= 1

2 3)

(18.18a)

ただしヤコビ行列A

i

=

@F

i

=@qは次のように与えられる

5

． A

i

=

0 B B B B B B B B B B B B B @

0 i

1 i

2 i

3

0 0 0

;u

1

u

i

+

i

1

2

D

i

1 i

2

u

1

;

i

1 ~

u

2 i

3

u

1

;

i

1

u

~

3 i

1 ~

0 0

;u

2

u

i

+

2 i

1

u

2

;

i

2

u

~

1

D

i

2 i

3

u

2

;

i

2

u

~

3 i

2 ~

0 0

;u

3

u

i

+

3 i

2 i

1

u

3

;

i

3

u

~

1 i

2

u

3

;

i

3 ~

u

2

D

i

3 i

3 ~

0 0

;Hu

i

+

2

u

i

1

H;

~

u

1

u

i

2

H;u

~

2

u

i

3

H;

~

u

3

u

i

u

i

0 0

;ku

i

1

k

i

2

k

i

3

k

0

u

i

0

;"u

i

1

"

i

2

"

i

3

"

0

u

i

(18.18b)

D

ij

=

u

i

+

ij

(1

;

~

)

u

j

~

=

;

1

2 = ~

uuu

2

=

2

である．なお次の関係がある． H

=

e

+

p

=

c

2 +

2 ~

=

e ;

2

このヤコビ行列A

i

は簡単とは言えない．もちろんこれを直接対角化することもできるが，ここでは非保存形方程式のものが簡単になることに着目しこれを利用して対角化を行う．式

(18.15)

に対応する非保存形方程式は次のようになる． @ @t

+

uuur

+

ruuu

= 0

(18.19a)

@ u u u @t

+

u u uru u u

+ 1

rp

= 0

(18.19b)

@p @t

+

uuurp

+

c

2

ruuu

= 0

(18.19c)

@k @t

+

uuurk

= 0

(18.19d)

@" @t

+

u u ur"

= 0

(18.19e)

この方程式系は

1

階の準線形微分方程式で，１つの空間次元の式はベクトル形で次のように書くことができる． @q @t

+

A

i

@q @x

i

= 0

(

i

= 1

2 3)

(18.20a)

5 Aの要素は，まずFの成分の式をqの成分で表し，それからqの成分で偏微分すれば容易に求めることができる．

(12)

ただし q

=

1

u

2

u

3

p k " 1 C C C C C C C C C C C C C A A

i

=

i

1 i

2 i

3

0 0 0

0

u

i

0

0 i

1

=

0 0

0

u

i

0 i

2

=

0 0

0

u

i

3

=

0 0

0 i

1

c

2 i

2

c

2 i

3

c

2

u

i

0 0

0

u

i

0

u

i

(18.20b)

式

(18.20)

は dq

=

@q @q dq

=

Ndq

(18.21)

であるから，式

(18.18)

の左からNを演算したものである．なお A

i

=

NA

i

N ;

1

であるが，q はqの成分の

1

次の同次式ではないから，上記のようにq

=

NqとdNq

= 0

は成立しない．保存形方程式を非保存形に変換する行列N

=

@q =@qとその逆行列N ;

1

は次のようになる

6

． N

=

1

0

0 0 0

0

;u

1

=

1

=

0

0 0 0

0

;u

2

=

0

1

=

0 0 0

0

;u

3

=

0

1

=

0 0

0 ~

u u u

2

=

2

;u

~

1

;u

~

2

;u

~

3 ~

0

;k=

0

0 0 1

=

0

;"=

0

0 0 0 1

= 1 C C C C C C C C C C C C C A

(18.22a)

N ;

1 =

1

0

0 0 0

u

1

0

0 0 0

u

2

0

0 0 0

u

3

0

0 0 0

u u u

2

=

2

u

1

u

2

u

3

1

=

~

0 0

k

0

"

0

0 (18.22b)

２つの行列A

i

とA

i

を比較すれば，A

i

が単純であることは一目瞭然である．ここではまずA

i

を対角化し，A

i

の対角化は次式によって行う． A

i

=

N ;

1

A

i

N

=

N ;

1

R

i

L

i

N

=

R

i

L

i

(18.23)

6 Nの要素は，まずq の成分の式をqの成分で表し，それからqの成分で偏微分すれば容易に求めることができる．

(13)

i

は行列A

i

の固有値

ik

からなる対角行列，L

i

はA

i

の左固有ベクトルl

ik

からなる行列，R

i

=

L ;

1 i

は右固有ベクトルr

ik

からなる行列である．ここでは次にA

i

の固有値

ik

と左固有ベクトルl

ik

の求め方について簡単に説明する

(

固有値，固有ベクトルに関してはは前章参照

)

．固有値

ik

は

i

の

7

次代数方程式 jA

i

;

i

Ij

=

u

i

;

i

1 i

2 i

3

0

u

i

;

i

0

0 i

1

=

0

u

i

;

i

0 i

2

=

0

u

i

;

i

3

=

0

0 i

1

c

2 i

2

c

2 i

3

c

2

u

i

;

i

0

u

i

;

i

0

u

i

;

i

= 0

の根で，

i

1 =

u

i

(5

重根

)

i

2 =

u

i

+

c

i

3 =

u

i

;cとなる．左固有ベクトルl

ik

を求める式は固有値

i

=

u

i

に対しては l

i

(

A

i

;

i

I

) = (

l

1 ik

l

2 ik

l

7 ik

)

0 i

1 i

2 i

3 0 0 0

0

0 i

1

=

0 0

0

0 i

2

=

0 0

0

0 i

3

=

0 0

0 i

1

c

2 i

2

c

2 i

3

c

2 0 0 0

0

0 0 0 0

0

0 0 0 0

= 0

固有値

i

=

u

i

cに対しては l

i

(

A

i

;

i

I

) = (

l

1 ik

l

2 ik

l

7 ik

)

0 B B B B B B B B B B B B B @ c

i

1 i

2 i

3

0 0 0

0

c

0

0 i

1

=

0 0

0

c

0 i

2

=

0 0

0

c

i

3

=

0 0

0 i

1

c

2 i

2

c

2 i

3

c

2

c

0 0

0

c

0

c 1 C C C C C C C C C C C C C A

= 0

となる．これらの式から

1

次独立の

7

個の左固有ベクトルを求めることができる．対角行列

i

の固有値と行列L

i

の固有ベクトルの順番はある任意性を持つが，ここではこの対角化によって得られる常微分方程式の系が，上から流速u

i

で伝播するエントロピー波，３つの，速度u

i

+

cで伝播する圧力波または流速u

i

で伝播する剪断波，速度u

i

;cで伝播する圧力波，流速u

i

で伝播するkの波，流速 u

i

で伝播する"の波の式になるように並べている．固有値と固有ベクトルは対応するところに配さなければならない．このときL

i

の対角要素は

0

でないある大きさを持つ．右固有ベクトルの行列R

i

は，

i

とL

i

の任意性のためL

i

の逆行列として求められる．行列A

i

の固有値も

i

である．固有値は波の速度を表すので，同一物理現象ではそれを記述する数式によらず同じになる．A

i

の左固有ベクトルはL

i

=

L

i

N，右固

(14)

有ベクトルはR

i

=

N ;

1

R

i

で，これらの行列は次のように表される．

i

=

i

0

0 0 0

0

u

i

+

i

1

c

0

0 0 0

0

u

i

+

i

2

c

0

0 0 0

0

u

i

+

i

3

c

0 0 0

0

u

i

;c

0 0

0

u

i

0

u

i

(18.24a)

L

i

=

1 0 0 0

;

1

=c

2 0 0

0 1 0 0

i

1

=c

0 0

0 0 1 0

i

2

=c

0 0

0 0 0 1

i

3

=c

0 0

0 i

1 i

2 i

3

;

1

=c

0 0

0 0 0 0

0 1 0

0 0 0 0

0 0 1

1 C C C C C C C C C C C C C A N

(18.24b)

R

i

=

N ;

1

1 i

1

=

2

c

i

2

=

2

c

i

3

=

2

c ;=

2

c

0 0

0 1

;

i

1

=

2

0

0 i

1

=

2 0 0

0

1

;

i

2

=

2

0 i

2

=

2 0 0

0

1

;

i

3

=

2 i

3

=

2 0 0

0 i

1

c=

2 i

2

c=

2 i

3

c=

2

;c=

2 0 0

0

0 1 0

0

0 0 1

(18.24c)

非保存形の式

(18.20)

と等価な常微分方程式系は，この式の左から左固有ベクトルL

i

を乗じることによって得られるが，これは保存形の式

(18.16)

の左から左固有ベクトルL

i

を乗じたものでもある． L

i

@q @t

+

@F

i

@x

i

=

L

i

@q @t

+

i

L

i

@q @x

i

= 0

(18.25)

この常微分方程式系は具体的に書けば次のようになる． ; @ @t

+

u

i

@ @x

i

;

1

c

2

; @ @t

+

u

i

@ @x

i

p

= 0

n @ @t

+(

u

i

+

ij

c

)

@ @x

i

o u

j

+

ij

c n @ @t

+(

u

i

+

c

)

@ @x

i

o p

= 0

(

j

= 1

2 3)

(18.26)

n @ @t

+(

u

i

;c

)

@ @x

i

o u

i

;

1

c n @ @t

+(

u

i

;c

)

@ @x

i

o p

= 0

; @ @t

+

u

i

@ @x

i

k

= 0

; @ @t

+

u

i

@ @x

i

"

= 0

これらの式は特性曲線dx

i

=dt

=

u

i

dx

i

=dt

=

u

i

+

c dx

i

=dt

=

u

i

; cに沿う内微分@=@t

+

u

i

@=@x

i

@=@t

+(

u

i

+

c

)

@=@x

i

または@=@t

+(

u

i

;c

)

@=@x

i

のいずれか１つを含むので常微分方程式で，上からdx

i

=dt

=

u

i

に沿って伝播するエントロピー波，dx

i

=dt

=

u

i

に沿う剪断波またはdx

i

=dt

=

u

i

+

cに沿う圧力波，dx

i

=dt

=

u

i

;c に沿う圧力波，dx

i

=dt

=

u

i

に沿う乱れの量kと"の波を表している．

(15)

18.2.2

一般曲線座標系の方程式

一般曲線座標系のオイラー方程式とこれに加えられるkと"の輸送方程式は次のようになる． @ @t J

+

@ @

i

JU

i

= 0

(18.27a)

@ @t Juuu

+

@ @

i

J

(

uuuU

i

+

r

i

p

) = 0

(18.27b)

@ @t Je

+

@ @

i

JHU

i

= 0

(18.27c)

@ @t Jk

+

@ @

i

JkU

i

= 0

(18.27d)

@ @t J"

+

@ @

i

J"U

i

= 0

(18.27e)

この方程式系から１つの空間次元の部分のみを取り出せば次式が得られる． @q

^

@t

+

@F

^

i

@

i

= 0

(

i

= 1

2 3)

(18.28a)

ただし

^

q

=

J 0 B B B B B B B B B B B B B @ u

1

u

2

u

3 ^

i

=

J 0 B B B B B B B B B B B B B @ U

i

u

1

U

i

+

i

1

p u

2

U

i

+

i

2

p u

3

U

i

+

i

3

p HU

i

kU

i

"U

i

(

i

= 1

2 3)

(18.28b)

流束F

^

i

はq

^

の成分だけでなく測度

ij

の関数であるが，時空間内の固定された点を考えれば，F

^

i

はq

^

の成分だけの

1

次の同次式で

Euler

の同次関係を用いれば dF

^

i

=

@F

^

i

@q

^

dq

^

= ^

A

i

dq

^

F

^

i

= ^

A

i

q

^

dA

^

i

q

^

= 0

(18.29)

となり，式

(18.28)

は次のよう線形化できる． @q

^

@t

+ ^

A

i

@q

^

@

i

= 0

(

i

= 1

2 3)

(18.30a)

^

i

=

@F

^

i

=@q

^

は次のように与えられる．

^

A

i

=

0 i

1 i

2 i

3

0

;u

1

U

i

+

i

1

2

D

^

1 i

2

u

1

;

~

i

1

u

2 i

3

u

1

;

~

i

1

u

3 ~

i

1

0

;u

2

U

i

+

i

2

2 i

1

u

2

;

~

i

2

u

1

D

^

2 i

3

u

2

;

~

i

2

u

3 ~

i

2

0

;u

3

U

i

+

i

3

2 i

1

u

3

;

~

i

3

u

1 i

2

u

3

;

~

i

3

u

2

D

^

3 ~

i

3

0

;HU

i

+

2

U

i

1

H;

~

U

i

u

1 i

2

H;U

~

i

u

2 i

3

H;

~

U

i

u

3

U

i

0

;kU

i

1

k

i

2

k

i

3

k

0

U

i

0

;"U

i

1

"

i

2

"

i

3

"

0

U

i

(18.30b)

(16)

^

D

j

=

U

i

+(1

;

~

)

ij

u

j

~

=

;

1

2 = ~

uuu

2

=

2 ij

=

@

i

=@x

j

である．デカルト座標系の非保存形方程式

(18.19)

は，uuur

=

U

i

@=@

i

r

=

r

i

@=@

i

であるから，容易に一般曲線座標系の式に書換えることができ次のようになる． @ @t

+

U

i

@ @

i

+

r

i

@ uuu @

i

= 0

(18.31a)

@ uuu @t

+

U

i

@ uuu @

i

+ 1

r

i

@p @

i

= 0

(18.31b)

@p @t

+

U

i

@p @

i

+

c

2

r

i

@ u u u @

i

= 0

(18.31c)

@k @t

+

U

i

@k @

i

= 0

(18.31d)

@" @t

+

U

i

@" @

i

= 0

(18.31e)

^

@t

+ ^

A

i

@q

^

@

i

= 0

(

i

= 1

2 3)

(18.32a)

ただし

^

q

=

q

=

0 B B B B B B B B B B B B B @ u v w p k " 1 C C C C C C C C C C C C C A A

^

i

=

0 B B B B B B B B B B B B B @ U

i

1 i

2 i

3

0 0 0

0

U

i

0

0 i

1

=

0 0

0

U

i

0 i

2

=

0 0

0

U

i

3

=

0 0

0 i

1

c

2 i

2

c

2 i

3

c

2

U

i

0 0

0

U

i

0

U

i

(18.32b)

非保存形方程式

(18.32)

は保存形方程式

(18.30)

の左から@q

^

=@q

^

=

N=Jを演算することによって得ることもできる．dq

^

= (

@q

^

=@q

^

)

dq

^

= (

N=J

)

dq

^

．式

(18.32)

は双曲型で，行列A

^

i

は対角化することができ，特性の理論により常微分方程式の系に変換できる．行列A

^

i

は多くの

0

要素を含むので容易に対角化でき，これをもとにA

^

i

も対角化できる．

^

A

i

=

N ;

1 ^

A

i

N

=

N ;

1 ^

R

i

^

i

L

^

i

N

= ^

R

i

^

i

L

^

i

(18.33)

行列A

^

i

の固有値と固有ベクトルはデカルト座標系の場合と同様の方法で求めることができる．その固有値は

^

i

1 =

U

i

(5

重根

)

，

^

i

2 =

U

i

+

p g

ii

c

^

i

3 =

U

i

; p g

ii

cとなる．g

ij

=

ik

jk

．行列A

^

i

の固有値の対角行列，左右の固有ベクトルの行列は次のようになる．

^

i

=

i

0

0 0 0

0

U

i

+

1 i

p g

11

c

0

0 0 0

0

U

i

+

2 i

p g

22

c

0

0 0 0

0

U

i

+

3 i

p g

33

c

0 0 0

0

U

i

; p g

ii

c

0 0

0

U

i

0

U

i

(18.34a)

(17)

^

L

i

=

1

0

;

1

=c

2 0 0

0 ii

1 i

1

2

;

2 i

2

1

1 i

1

3

;

3 i

3

1

1 i

p g

11

=c

0 0

0

2 i

2

1

;

1 i

1

2 ii

2 i

2

3

;

3 i

3

2

2 i

p g

22

=c

0 0

0

3 i

3

1

;

1 i

1

3

3 i

3

2

;

2 i

2

3 ii

3 i

p g

33

=c

0 0

0 i

1 i

2 i

3

; p g

ii

=c

0 0

0

0 1 0

0

0 0 1

1 C C C C C C C C C C C C C A N

(18.34b)

^

R

i

=

N ;

1

0 B B B B B B B B B B B B B B B B B B B B B B @

1

1 i

2

p g

11

c

2 i

2

p g

22

c

3 i

2

p g

33

c ;

2

p g

ii

c

0 0

0

D

1

3

2

2 i

;

1 i

2

g

ii

3

2

3 i

;

1 i

3

g

ii

i

1

2

g

ii

0 0

0

3

2

1 i

;

1 i

2 i

1

g

ii

D

2

3

2

3 i

;

1 i

2 i

3

g

ii

i

2

g

ii

0 0

0

3

2

1 i

;

1 i

3 i

1

g

ii

3

2

2 i

;

1 i

3 i

2

g

ii

D

3 i

3

2

g

ii

0 0

0

1 i

c

2

p g

11

2 i

c

2

p g

22

3 i

c

2

p g

33

; c

2

p g

ii

0 0

0

0 1 0

0

0 0 1

1 C C C C C C C C C C C C C C C C C C C C C C A

(18.34c)

D

j

=

3

2 ji

;

1 ij

2

g

ii

+(1

;

ji

) 1

ii

(

j

= 1

2 3)

式

(18.28)

と等価な常微分方程式系は次のようになる．

^

L

i

J @q

^

@t

+

@F

^

i

@

i

= ^

L

i

@q

^

@t

+ ^

i

L

^

i

@q

^

@

i

= 0

(18.35)

あるいは @ @t

+

U

i

@ @

i

;

1

c

2

@ @t

+

U

i

@ @

i

p

= 0

n @ @t

+

; U

i

+

p g

ii

c @ @

i

o U

i

+

p g

ii

c n @ @t

+

; U

i

+

p g

ii

c @ @

i

o p

= 0

(

j

=

i

)

ii

@ @t

+

U

i

@ @

i

u

j

;

ij

@ @t

+

U

i

@ @

i

u

i

= 0

(

j6

=

i

)

(18.36)

n @ @t

+

; U

i

; p g

ii

c @ @

i

o U

i

; p g

ii

c n @ @t ; ; U

i

; p g

ii

c @ @

i

o p

= 0

@ @t

+

U

i

@ @

i

k

= 0

@ @t

+

U

i

@ @

i

"

= 0

これらの式は空間内の特性曲線d

i

=dt

=

U

i

d

i

=dt

=

U

i

+

p g

ii

c d

i

=dt

=

U

i

; p g

ii

cに沿う内微分 @=@t

+

U

i

@=@

i

@=@t

+(

U

i

+

p g

ii

c

)

@=@

i

@=@t

+(

U

i

; p g

ii

c

)

@=@

i

のいずれか１つのみを含む常微分方程式で，写像空間の特性曲線に沿って伝播するエントロピー波，圧力波，剪断波，kの波，"の波のいずれかを表している．

(18)

18.2.3

一般曲線座標系の質量流束の方程式

一般曲線座標系の質量流束のオイラー方程式とこれに加えられるkと"の方程式は次のようになる． @ @t J

+

@ @

i

JU

i

= 0

(18.37a)

@ @t JU U U

+

@ @

i

J ; U U UU

i

+(

r

i

r

)

p

= 0

(18.37b)

@ @t Je

+

@ @

i

JHU

i

= 0

(18.37c)

@ @t Jk

+

@ @

i

JkU

i

= 0

(18.37d)

@ @t J"

+

@ @

i

J"U

i

= 0

(18.37e)

~

@t

+

@F

~

i

@

i

= 0

(

i

= 1

2 3)

(18.38a)

ただし

~

q

=

1

U

2

U

3 ~

i

=

i

U

1

U

i

+

g

i

1

p U

2

U

i

+

g

i

2

p U

3

U

i

+

g

i

3

p HU

i

kU

i

"U

i

(

i

= 1

2 3)

(18.38b)

流束F

~

i

はq

~

の成分の

1

次の同次式であるから，オイラーの同次関係を用いれば dF

~

i

=

@F

~

i

@q

~

dq

~

= ~

A

i

dq

~

F

~

i

= ~

A

i

q

~

dA

~

i

q

~

= 0

(18.39)

となり，式

(18.38)

は次のよう線形化できる． @q

~

@t

+ ~

A

i

@q

~

@

i

= 0

(

i

= 1

2 3)

(18.40a)

~

i

=

@F

~

i

=@q

~

は

~

A

i

=

0

1 i

2 i

3 i

0 0 0

;U

1

U

i

+

g

i

1

2

D

~

i

1

2 i

U

1

;g

~

i

1

2

3 i

U

1

;

~

g

i

1

3

g

~

i

1 0 0

;U

2

U

i

+

g

i

2

1 i

U

2

;

~

g

i

2

1

D

~

i

2

3 i

U

2

;

~

g

i

2

3

g

~

i

2 0 0

;U

3

U

i

+

g

i

3

2

1 i

U

3

;

~

g

i

3

1

2 i

U

3

;

~

g

i

3

2

D

~

i

3

g

~

i

3 0 0

;HU

i

+

2

U

i

1 i

H;

~

1

U

i

2 i

H;

~

2

U

i

3 i

H;

~

3

U

i

U

i

0 0

;kU

i

1 i

k

2 i

k

3 i

k

0

U

i

0

;"U

i

1 i

"

2 i

"

3 i

"

0

U

i

(18.40b)

(19)

~

D

ij

=

U

i

+

ji

U

j

;g

~

ji

j

i

= (

@x

k

=@

i

)

u

k

=

g

ij

U

j

g

ij

= (

@x

k

=@

i

)(

@x

k

=@

j

)

．この場合の非保存形方程式は次のようになる． @ @t

+

U

i

@ @

i

+

@U

i

@

i

= 0

(18.41a)

@ UUU @t

+

U

i

@ UUU @

i

+ 1

r

i

r @p @

i

= 0

(18.41b)

@p @t

+

U

i

@p @

i

+

c

2

@U

i

@

i

= 0

(18.41c)

@k @t

+

U

i

@k @

i

= 0

(18.41d)

@" @t

+

U

i

@" @

i

= 0

(18.41e)

~

@t

+ ~

A

i

@q

~

@

i

= 0

(18.42a)

ただし

~

q

=

1

U

2

U

3

p k " 1 C C C C C C C C C C C C C A A

~

i

=

i

1 i

2 i

3 i

0 0 0

0

U

i

0

g

i

1

=

0 0

0

U

i

0

g

i

2

=

0 0

0

U

i

g

i

3

=

0 0

0

1 i

c

2

2 i

c

2

3 i

c

2

U

i

0 0

0

U

i

0

U

i

(18.42b)

この非保存形方程式

(18.42)

は，保存形方程式

(18.40)

の左から行列@q

~

=@q

~

= ~

N=Jを演算することによっても得られる．保存形から非保存形への変換の行列N

~

=

J@q

~

=@q

~

とその逆行列N

~

;

1

は次のようになる．

~

N

=

1

0

0 0 0

0

;U

1

=

1

=

0

0 0 0

0

;U

2

=

0

1

=

0 0 0

0

;U

3

=

0

1

=

0 0

0

2

;

~

1

;

~

2

;

~

3 ~

0

;k=

0

0 0 1

=

0

;"=

0

0 0 0 1

= 1 C C C C C C C C C C C C C A

(18.43a)

~

N ;

1 =

1

0

0 0 0 0

U

1

0

0 0 0 0

U

2

0

0 0 0 0

U

3

0

0 0 0 0

uuu

2

=

2

1

2

3

1

=

~

0 0

k

0

"

0

0 0 0

(18.43b)

(20)

なおdq

~

= (

@q

~

=@q

~

)

dq

~

= ( ~

N=J

)

dq

~

．行列A

~

i

は簡単で容易に対角化でき，これよりA

~

i

も対角化できる．

~

A

i

= ~

N ;

1 ~

A

i

N

~

= ~

N ;

1 ~

R

i

~

i

L

~

i

N

~

= ~

R

i

~

i

L

~

i

(18.44)

行列A

~

i

の固有値の対角行列

~

i

，その左固有ベクトルの行列L

~

i

，右固有ベクトルの行列R

~

i

= ~

L ;

1 i

は次のようになる．

~

i

= ^

i

(18.45a)

~

L

i

=

1

0

;

1

=c

2 0 0

0

1

;

2 i

g

12

=g

22

;

3 i

g

13

=g

33

1 i

p g

11

=c

0 0

0

;

1 i

g

21

=g

11

1

;

3 i

g

23

=g

33

2 i

p g

22

=c

0 0

0

;

1 i

g

31

=g

11

;

2 i

g

32

=g

22

1

3 i

p g

33

=c

0 0

0

1 i

2 i

3 i

; p g

ii

=c

0 0

0

0 1 0

0

0 0 1

~

N

(18.45b)

~

R

i

= ~

N ;

1

1 i

=

2

p g

11

c

2 i

=

2

p g

22

c

3 i

=

2

p g

33

c ;=

2

p g

ii

c

0 0

0

1

;

1 i

=

2

2 i

g

12

=

2

g

22

3 i

g

13

=

2

g

33

g

1 i

=

2

g

ii

0 0

0

1 i

g

21

=

2

g

11

1

;

2 i

=

2

3 i

g

23

=

2

g

33

g

2 i

=

2

g

ii

0 0

0

1 i

g

31

=

2

g

11

2 i

g

32

=

2

g

22

1

;

3 i

=

2

g

3 i

=

2

g

ii

0 0

0

1 i

c=

2

p g

11

2 i

c=

2

p g

22

3 i

c=

2

p g

33

;c=

2

p g

ii

0 0

0

0 1 0

0

0 0 1

(18.45c)

式

(18.38)

と等価な常微分方程式系は次のようになる．

~

L

i

J @q

~

@t

+

@F

~

i

@

i

= ~

L

i

@q

~

@t

+ ~

i

L

~

i

@q

~

@

i

= 0

(18.46)

なおこの式は，常微分方程式系に書き替えれば当然のことながら式

(18.36)

になり，写像空間の特性曲線に沿って伝播するエントロピー波，圧力波，剪断波，kと"の波の式である．