オイラー方程式の数値解

(1)

オイラー方程式の数値解

Numerical solution of Euler equation

物理学専攻　宮崎巧也

Department of physics Miyazaki Takuya

完全流体については，オイラー方程式で記述される．この方程式の解が有限時間内に無くなるかどうかという疑問が数学者，物理学者の間で現在も存在する．つまり，物理量の時間発展を追跡した場合，この物理量が有限時間内に発散するのではないか，という疑問が未に残っている．しかし，この方程式は強い非線形性を持っていて，一般にその解析解を求めるのは困難である．そのため，解を求めるには数値計算を用いる．この研究を始めた動機について述べておく。４次元のナヴィエ＝ストークス方程式（以後簡単のために

NS

方程式と略記する）のシミュレーション

[1],[2]

、から分かっていることは、（１）３次元と同様にエネルギーカスケードが起こる、（２）波数空間でのエネルギー輸送は３次元より速いということであった。３次元

NS

方程式では渦糸同士の衝突とつなぎ変えがエネルギー散逸を増加させることが確かめられている

[3]

。４次元で散逸がより大きいのは衝突の頻度がより大きいことを意味するが，特異構造のトポロジーはどうであろうか？粘性のないオイラー方程式で、衝突後の解の継続的な存在が困難視されることから、特異構造の衝突の効果はもっと重要になるであろう。そのため３次元と４次元のオイラー方程式の解の振る舞いと特異構造のトポロジーを比較検討することは、オイラー方程式の解の発散の問題の解明に役立つ。オイラー方程式で記述される非粘性流体では滑らかな速度場から出発しても、方程式の移流項が大きなスケールの揺らぎを小さなスケールに輸送するので、特徴的な長さスケール

l(t)

は時間と共に減少するであろう．

l(t)

は指数的にしか減少しないのか、

それとも有限時間でゼロになるのか。オイラー方程式ではエネルギーは保存されるから、

代表的な速度の大きさ

u

0 は一定のままであるが、その空間微分は

u

0

/l(t)

のようにスケールされるので、任意の次元で定義される速度場の

2-form

である

Ω

ij

= ∂u

j

∂x

i

− ∂u

i

∂x

j

の振幅が時間と共にどのように増大するかが研究対象になる

[4]

。通常は

Ω

ij の２乗である

Ω

2

= Σ

ij

Ω

²_ij

で定義されたエンストロフィーの時間発展が調べられる。

Morf

^等

[5]

^は

Taylor-Green

^流れの初期場に対して，

Ω

2

(t)

を時間のベキ級数で求めて、その発散の具合を調べた。

Kerr [6]

は互いに接近する反対向きの渦度を持つ２本の渦糸系が発散に導くことを示し、

Pelz

等

[7]

は対称性の高い

Taylor-Green

流れを初期場として発散を確かめた。

Ohkitani

等

[8]

は

Burgers

渦を拡張したモデル方程式の数値解析により有限時間の発散を確かめた。もし発

散が起こるとすれば、その原因はなにか。渦がストレインによって引き伸ばされて、渦度が増大するのは一つの原因である。あるいは反対の極性を持つ特異構造が衝突を起こし、

速度の空間的勾配が増大するメカニズムも原因と考えられる。もし粘性があれば、渦のつなぎかえが起こるが、粘性がなければそのようなことが起こらずに破局を迎えるだろう。

ここでオイラー方程式と

NS

方程式での特異構造のトポロジーの違いについて述べておく必要がある。３次元では

Kelvin-Helmholtz

不安定によりシート状の渦が巻き上がる。粘性がなければシートはますます巻き込まれ、大きなスケールで眺めれば、それは線状に見えるであろう。粘性があれば、シート間の速度勾配が大きくなり、粘性の効果でそれらは均されて線になる。ある程度のスケールで眺めれば、粘性があろうとなかろうと構造は同じ線状であると考えられ、我々の行った３次元でのオイラー方程式と

NS

方程式の渦度の可

(2)

視化もその描像を支持している。微小スケールでの構造の同定は

NS

方程式における方が易しい。なぜならオイラー方程式では最小スケールがメッシュサイズよりも小さくなると、

微細構造に悪影響が現れ、その同定が難しいからである。

NS

方程式では乱流が減衰してしまわない限り微細構造は安定に存在し続ける。４章で述べる特異構造のトポロジーの同定は、粘性の小さな高レイノルズ数の

NS

方程式のものとオイラー方程式のそれとは同じであると仮定し、

NS

方程式のシミュレーションの結果に基づいている。

４章では４次元乱流の特異構造のトポロジーについて議論している．特異構造の繋ぎ換えが時空間で起こる条件は、空間次元

d

と特異構造の次元

d

s とすれば、

2(d

s

+ 1) ≥ d + 1

である。発達した３次元古典乱流では特異構造は渦糸

(d

s

= 1)

であるので、この条件が辛うじて満たされる。もし４次元で

d

s

= 1

なら、繋ぎ換えは起こらず、流れは乱流にならない。４次元の特異構造の次元はいくらであろうか。本研究では、

(1)

渦度の一般化である

2-form Ω

ij

= ∂u

j

/∂x

i

− ∂u

i

/∂x

jの可視化、

(2)

_ij

Ω

²_ij

, Ω

²12の一般化次元、

(3)

ストレインテンソルの固有値の分布、を調べることにより、４次元の特異構造の次元は

2

^であることを確かめた。まず

128

⁴のシミュレーションで、

ρ =

_ij

Ω

²_ij の可視化を行った。第４の座標

x

4を固定した３次元空間での構造は、３次元乱流の渦糸構造とよく似て、ライン構造が見られる。

x

4を変えても、構造が徐々に変わるだけなので、４次元の特異構造はシート状であることを示している。次にマルチフラクタル解析を用いて

ρ

の一般化次元

D

qとその１個の成分

Ω

²12の一般化次元

D ˜

qを調べた。レイノルズ数が低いので、右図のように大きな

q

^でほぼ

2.7

^{程度でシートの次元}

2

より少し大きい。同じレイノルズ数の３次元乱流と比べて、右図のように４次元乱流の一般化次元は３次元のそれプラス１であることが確かめられた。発達した３次元乱流の特異構造は線であるので、４次元では面であることと一致している。ストレインテンソルの固有値分布は２伸張１圧縮を示しており、１伸長１圧縮のライン構造を持つ３次元乱流と比べることにより、４次元はシート構造であることを示唆する。以上，４次元乱流の特異構造はシートであり、特異構造のつなぎかえは起こり、乱流は実現されることが確認できた。

５章ではオイラー方程式の数値解について議論している．オイラー方程式において、有限時間で特徴的なスケール

l(t)

がゼロになり、エンストロフィーが発散するかどうかは流体力学で最も重要な問題の一つである。

DNS

でオイラー方程式を扱う際の問題点は、特徴的な長さ

l(t)

がメッシュサイズより小さくなってしまうことである。それに対処する通常の方法では、エンストロフィーを時間のベキ級数で展開してその発散の具合を調べる

[1]

。本報告ではエンストロフィーをメッシュサイズを取り入れたスケーリング関数で表し、発散の問題を調べた。本研究では、エネルギースペクトルの

n

次モーメント

Ω

n

=

E(k)k

ⁿ

dk

のうち、

Ω

2 に注目する。典型的な速度を

u

0、スケールを

l(t)

、メッシュサイズを

a

とすれば、

Ω

2

(t) = u

²0

/l(t)

²

f (a/l(t))

と表される。

u

²0

= Ω

0を代入すれば、規格化されたエンストロフィー

R(t) = Ω

2

/Ω

0は

R(t) = 1 l(t)

²

f

a

l(t)

(1)

のようにスケーリング関数

f

を用いて、

l(t), a

だけで表される。

a → 0

の極限では、

R(t) ∼

l(t)

⁻²のような振る舞いが期待される。ほぼ同じ初期速度場から出発した種々のメッシュ数

N = 32

³

, 64

³

, 128

³

, 256

³

, 512

³ での

R(t)

をプロットすると，カットオフの影響で飽和していくが，

N → ∞

での振舞いは，それらの曲線群の包絡線であり、下に凸である。明らかに包絡線は指数関数より速いので、

R(t)

は指数関数より速く増大することは確かである。

l(t)

の振る舞いの定量的な評価を得るためには、

(1)

を

R(t)/a

²

= (a/l(t))

²

f (a/l(t)) = g(a/l(t))

と書き換えると、

R(t)/a

²は唯一の変数

a/l(t)

の関数値である。

l(t) = l

0

(t

_∗

− t)

⁻ⁿと置くことにより、３個のデータセットから

n, t

_∗を決定することが出来る。その結果を図２に掲げる。どのようなデータセットを用いても、ある時間の後、

n, t

_∗共に一定の値になるはず

(3)

図

1:

３次元に於ける各時刻の

n, t

_∗

である。図１を見ると，数値計算で得られた

n, t

_∗は

k

⁰

u

⁰

t = 1.0 ∼ 2.0

の区間で一定である．この事は導入したスケーリング関数がうまく機能している証拠である．

結論をまとめると，４次元乱流の特異構造がシート状であること，それに伴い，４次元乱流では特異構造の衝突頻度が３次元よりも高いことが分かった．また，オイラー方程式でエンストロフィー

Ω

2の時間発展を追跡したところ４次元の方が速く増大する傾向を観測できた．２次元は点状，３次元は線状，４次元は面状の特異構造を持つことから帰納的に考えれば，特異構造の次元

d

sと乱流の次元

d

の間に

d

s

= d − 2

という関係が成り立つ事を予測できるが５次元では満たされているであろうか？それに伴い，次元を増やして行けば

Ω

2は，より速く増大することが予想できるが実際はどうであろうか？本研究では３次元，４次元に注目したが偶数次元は奇数次元と違い，非粘性不変量が無限個存在するため，系に何らかの束縛を与えている可能性もあり，今後５次元乱流を詳しく調べ３次元乱流と比較することでより平等な比較が可能となるであろう．

また

Ω

²のメッシュ幅依存性を考慮したモデルにより，メッシュ幅

a

^→

0

^{の極限での振る} 舞いを予測した．先行の研究で

δ(t)

がメッシュ幅に到達するまでの時刻（本研究で言う所

の

k

0

u

0

t = 0.63

）までしか扱えていなかった．しかし，本研究ではカットオフの影響（飽和

の具合）を取り入れることで

δ(t)

がメッシュ幅に到達した後の時刻の

n,t

_∗を得ることが出来た．

k

0

u

0

t = 1.0

と

k

0

u

0

t = 2.0

の間では本研究で導出したスケーリング則が成り立っており，先行研究と本研究は相補的な関係にある．

t

_∗については４次元で

k

0

u

0

t

_∗

= 1.9

であり，３次元の

k

0

u

0

t

_∗

= 2.8

^{と比べると発散時刻}

t

_∗が早い事がわかった．これは４次元の方が３次元に比べ特異構造が接近する頻度が高い事を示しておりトポロジーで行った考察と矛盾しない．また

n

については４次元では

n = 0.58

と読み取れるが３次元の

n = 0.62

に近い値となった．今後の発展のために改善すべき点を述べると，本研究で用いたメッシュ数では

n

，

t

_∗は時間的に僅かに揺らいでおり，極小での

n, t

_∗の値を信頼しても良いのか？

という疑問が残る．この揺らぎの原因がオイラー方程式の性質なのか，数値計算の精度のせいなのか，現段階では分からない．そのため今後，単純な１つのベキで記述できない可能性も視野に入れ，更に大規模な計算で同様の解析を行なう必要がある．４次元に関しては３次元と同じサイズで解析を行い両者を比較する必要がある．また，本研究に限った話では無いが

DNS

で非粘性不変量が保存している時間は初期からほんの僅かな時間だけで

(4)

あるため，その影響を評価するには更なる一様等方性が必要となり，その意味でも大規模計算は欠かせない．

参考文献

オイラー方程式の数値解