筑波大学大学院博士課程
システム情報工学研究科修士論文

Wavelet Taylor Galerkin 法による１次元 Euler 方程式の数値計算松浦大志

修士（工学）（構造エネルギー工学専攻）

指導教員 松島

亘志

２０１９年

３月

概要
Euler 方程式は圧縮性非粘性流れの基礎方程式であり、圧縮性流れのシミュレーション手法のベンチマークとしてよく用いられる。Euler 方程式は衝撃波解を持つため、衝撃波付近を精度よく再現できることが、計算手法に求められる性質である。本研究では Euler 方程式のための新しい計算手法として、Wavelet Galerkin 法およびその改良版である Wavelet Taylor Galerkin法に着目した。本論文では、これらの手法の特性を見出し、既存の代表的な方法に対する優位性を示すことを目的としている。そこで、1 次元問題に焦点を当て、線形の移流方程式と、
Euler 方程式の数値シミュレーションを行い、特定の条件下で、既存の差分法の結果に比べて良好な結果が得られることを確認した。

目次
第 1 章 はじめに ................................................................................................................................. 1
第 2 章 Wavelet Galerkin 法の概要 .................................................................................................... 3
2.1
Wavelet Galerkin 法の概要と既往の研究 ............................................................................ 3
2.2
本研究で用いるスキームの導出 ......................................................................................... 32.3
線形移流方程式による検証と課題点 ................................................................................. 6
第 3 章 Wavelet Taylor Galerkin 法の概要とスキームの検証 ......................................................... 9
3.1
Taylor Galerkin 法の概要....................................................................................................... 9
3.2
Wavelet Taylor Galerkin 法による移流方程式の計算結果 ............................................... 10
3.3
Wavelet Taylor Galerkin 法に対する安定性解析 ................................................................11
3.4
Godunov の定理とスキームの非線形化 ........................................................................... 13
第 4 章 Wavelet Taylor Galerkin 法による Euler 方程式の数値シミュレーション .................... 15
4.1
提案手法による 1 次元 Euler 方程式の離散化 ................................................................ 15
4.2
計算モデルの概要 ............................................................................................................... 164.3
計算結果 ............................................................................................................................... 17
4.4
陽的人工粘性項の改良とその効果 ................................................................................... 18
4.5
人工粘性項の経験的なパラメータに関する検討 ........................................................... 19
第 5 章 まとめ ................................................................................................................................... 21
5.1
得られた知見 ....................................................................................................................... 21
5.2
今後の課題 ........................................................................................................................... 21
謝辞
22
参考文献 ............................................................................................................................................... 23
付録 A Wavelet Taylor Galerkin 法が 2 次以上の精度を持つことの検証 ................................... 24
付録 B Lax Wavelet Galerkin 法 ....................................................................................................... 26
付録 C 多段階 Wavelet Galerkin 法 ................................................................................................. 27
付録 D 人工粘性項のパラメータに関する追加検討 ................................................................... 28

図表目次
図 2-1

coiflet スケーリング関数 : 𝜙 ................................................................................................................. 6

図 2-2

基底関数 coiflet による展開係数 : 𝑎𝑙

表 2-1

移流方程式計算条件 .............................................................................................................................. 7

図 2-3

Wavelet Galerkin 法によるサインカーブの移流計算(coif1) ............................................................... 7

図 2-4

Wavelet Galerkin 法による矩形波の移流計算 ...................................................................................... 7

図 3-1

基底関数 coiflet による 2 階の展開係数 .............................................................................................. 9

図 3-2

Wavelet Taylor Galerkin 法による矩形波の移流計算......................................................................... 10

図 3-3

基底関数の種類と格子点における誤差の二乗和 ..............................................................................11

図 3-4

Wavelet Taylor Galerkin 法 (coiflet 基底)の安定性解析 ..................................................................... 12

図 3-5

安定性解析の検証(coif1) ...................................................................................................................... 13

図 3-6

陽的人工粘性の付加 ............................................................................................................................ 14

図 4-1

1 次元衝撃波管 ..................................................................................................................................... 17

表 4-1

初期条件 ................................................................................................................................................ 17

図 4-2

1 次元衝撃波管問題の理論解 (𝑡 = 0.1) ............................................................................................ 17

表 4-2

計算条件 ................................................................................................................................................ 18

図 4-3

1 次元衝撃波管問題の計算結果 .......................................................................................................... 18

図 4-4

人工粘性項の計算方法による比較 ..................................................................................................... 19

図 4-5

各手法の精度の比較 ............................................................................................................................ 19

図 B-1

Lax Wavelet Galerkin 法による移流方程式の数値シミュレーション ............................................ 26

図 C-1

3 段階 Wavelet Galerkin 法による移流方程式の数値シミュレーション ....................................... 27

図 D-1

誤差の 𝜀 依存性（Δ𝑡/Δ𝑥 = 0.01, 0.005） ....................................................................................... 28(1)................................................................................................ 6

第1章

はじめに

液体や気体の流れを解析する分野を流体力学といい、その中でも特に流体の密度が変化す
る性質、すなわち圧縮性を考慮に入れて流れを解析する分野を圧縮性流体力学という。エンジンおよび、配管内部の高速流れや、航空機など流体中を高速で運動する輸送機器周りの流れは、圧縮性を考慮して解析する必要がある。また、流れが低速であっても、流体中の圧力波の伝播が重要な問題では、流体の圧縮性を考慮する必要がある。これらの流れ場を解析、予測することは各種機器の合理的な設計を可能とし、産業的に有益である。流れ場の解析には、大きく分けて理論解析、実験、数値シミュレーションの 3 つの観点からのアプローチが存在し、それぞれが互いを補い合って流体力学の枠組みを形成している。数値シミュレーシ
ョンは、理論的に解くことが困難な流れ場や、実験環境を整備するのが困難な流れ場の理解のためのツールとして、計算機の性能とともに発展してきた。
流体力学における数値シミュレーションの分野は CFD(Computational Fluid Dynamics : 数値
流体力学)ともよばれ、CFD スキームには堅牢性、精度、汎用性、簡便さといった性能が要求される。しかし、これらの要求すべてを同時に満たすことは困難であり、現在においても、
既存手法の改良や、新たな手法の提案がなされている。[1][2]
圧縮性 CFD において、解くべき基礎方程式は一般に、圧縮性を考慮した Navier - Stokes 方程式である。この方程式は、質量、運動量、エネルギーそれぞれの保存則から成る非線形偏微分方程式系である。本研究が対象とする Euler 方程式は圧縮性 Navier - Stokes 方程式から、粘性項と熱流束項を取り除いた方程式である。多くの圧縮性 CFD は Euler 方程式を数学的な基礎としており、新たなスキームの開発において、Euler 方程式の数値シミュレーションを実施することでスキームを検証することは不可欠なステップである。
Euler 方程式は、非線形方程式であるため、その解は衝撃波のような不連続面を持ちうることが知られており、そのような不連続面の近傍における精度と安定性を確保することが、圧縮性 CFD スキームに特有の条件である。この要求を満たすことができる新しい手法として、
本研究が注目したのが Wavelet Galerkin 法[3]である。Wavelet Galerkin 法は、Galerkin 法の一種であり、その数学的な基礎は有限要素法と同じく重み付き残差法である。Wavelet Galerkin 法
は、基底関数に、有界な台を持つ正規直交 Wavelet 基底を用いることにより、衝撃波のような不連続面を形成する非線形偏微分方程式の解を精度よく計算できると考えられる。Wavelet
Galerkin 法において、Wavelet 基底の正規直交性は、有限要素法のように連立一次方程式を解かなくて済むという利点をもたらす。この場合、時間方向の離散化は前進差分を用いると、離散化式が非常にシンプルになるが、これは精度や安定性の観点から適切ではない。そこで
時間方向の近似をより高次精度にするために、Taylor Galerkin 法[4]の考え方を導入した
Wavelet Taylor Galerkin 法[5][6]が存在する。しかし、Wavelet Taylor Galerkin 法の多くのスキー
ムでは、基底関数の性質に基づき、空間方向の分割を規則的な格子幅で行っている。一方、山口ら[7]は、基底関数に対して変数変換を施すことにより任意の格子幅での分割を可能としたが、時間方向の離散化については前進差分を用いている。また、Wavelet Galerkin 法により、
Euler 方程式を解いた例としては、浜井ら[8]の検討があげられるが、人工粘性のみで安定化を行っており、時間方向の離散化には前進差分を用いている。
そこで本研究では、山口らの方法を参考に、Taylor Galerkin 法の考え方を導入したスキームを開発し、Euler 方程式の数値シミュレーションを行うこととした。本論文の構成は以下のようである。まず、第 2 章で山口らの方法を参考に Wavelet Galerkin 法のスキームを導出する。
第 3 章では、山口らの方法に Taylor Galerkin 法の考え方を導入したスキームを用いて、線形1

移流方程式の数値シミュレーションを行い、提案手法の妥当性、信頼性を検証する。第 4 章
で提案手法を 1 次元圧縮 Euler 方程式に適用し、既存の手法と定量的に比較することで、本手法が有効な計算条件を見出し、第 5 章で全体をまとめる。

2

第2章
2.1

Wavelet Galerkin 法の概要
Wavelet Galerkin 法の概要と既往の研究

微分方程式の数値解法は差分法、有限要素法、スペクトル法[9]など様々ある。中でも、有限要素法やスペクトル法の基礎となっているのは重み付き残差法という考え方である。
Wavelet Galerkin 法は、重み付き残差法における基底関数に、Wavelet 解析で使用される基底関数を用いて離散化式を導出する方法である。
Wavelet 解析[10]は、Fourier 解析と同様に、信号データの周波数特性を調べる方法である。
Fourier 解析が信号を三角関数の重ね合わせで表現するのに対し、Wavelet 解析では、有界な台を持つ関数を使用する。それにより、Fourier 解析よりも周波数特性の時間変化や空間変化を捉えやすいという特徴がある。このような特徴により、Wavelet 解析は特異点の検出に優れている。そのため、Wavelet 解析で用いられる基底関数を用いたスキームは、衝撃波のような不連続面が形成される非線形現象を精度よく計算できると考えられている。
Wavelet Galerkin 法は、1990 年に Liandrat[3]などにより提案された。その後、Kumar らによって、Taylor Galerkin 法[4]の考え方を導入した、Wavelet Taylor Galerkin 法[5][6]が開発され、移流方程式をはじめとする双曲型方程式や、拡散方程式をはじめとする放物型方程式の数値
シミュレーションが行われている。これらは、次節で述べる基底関数のトゥースケール関係に基づき、規則的な格子幅で離散化を行っている。一方、山口ら[7]は、基底関数に変数変換を施すことにより、任意の格子幅で離散化することが可能な Wavelet Galerkin 法のアルゴリズ
ムを提案した。しかし、山口らの検討では、時間微分項の離散化には Euler 前進差分を用いている。そこで本研究では、山口らの方法に、Taylor Galerkin 法の考え方を導入したスキームを
用いる。

2.2

本研究で用いるスキームの導出

本節では、前節で述べた山口らのスキームの導出過程を紹介する。前述のように、Wavelet
Galerkin 法は、Galerkin 法の基底関数に、Wavelet 解析で用いられる基底関数を利用する手法であるため、基底関数の性質が重要となる。まず、次に述べる 3 つの性質を持った、有界な台を持つ関数 𝜙(𝑥) をスケーリング関数と定義する。[10]ここで、台とは、その区間以外で、関数が常にゼロとなることである。
条件 1

𝜙𝑗,𝑘 (𝑥) = 2−𝑗/2 𝜙(2−𝑗 𝑥 − 𝑘)
とおくとき、𝜙𝑗,𝑘 (𝑥)が正規直交系をなす。つまり、(1)∞〈𝜙𝑗,𝑘 (𝑥), 𝜙𝑗,𝑙 (𝑥)〉 ≡ ∫ 𝜙𝑗,𝑘 (𝑥)𝜙𝑗,𝑙 (𝑥)d𝑥 = 𝛿𝑘𝑙

(2)

−∞

条件 2
𝜙0,0 (𝑥) を、𝜙−1,𝑘 (𝑥) の線形結合で表せる。つまり、以下を満たす数列 𝑝𝑘 が存在する

𝜙0,0 (𝑥) = √2 ∑ 𝑝𝑘 𝜙(2𝑥 − 𝑘)𝑘

条件 3

𝜙 のフーリエ変換 𝜙̂(𝑤)が 𝑤 = 0 の近くで連続かつ 𝜙̂(0) = 1

3(3)

ここで、𝑗 と 𝑘 は任意の整数、𝛿𝑘𝑙 は Kronecker のデルタを表す。また、条件 2 のような関係を、トゥースケール関係といい、𝑝𝑘 をトゥースケール数列とよぶ。
つぎに、𝜙𝑗,𝑘 (𝑥) は正規直交関数系であるので、任意の関数𝑓(𝑥)を、𝜙𝑗,𝑘 (𝑥) で展開すると、
𝑓𝑗 = ∑ 𝑠𝑗,𝑘 𝜙𝑗,𝑘 (𝑥)

(4)

𝑘∈ℤ

のようにかける。式(4)では、添え字 𝑗 が残る。𝑗 は基底となるスケーリング関数の横幅を決める数であり、関数 𝑓(𝑥) の近似の解像度をあらわす。式(4)は、𝑗 → −∞ で 𝑓(𝑥) に収束する。𝑠𝑗,𝑘 は近似係数とよばれ、スケーリング関数の直交性より
𝑠𝑗,𝑘 = 〈𝑓, 𝜙𝑗,𝑘 (𝑥)〉(5)で求められる。このような、スケーリング関数のみを使う方法を、山口らの方法に倣い、𝜙 ス
キームとよぶ。今回の検討では、すべてこのスキームを使う。つぎに、𝑘 が整数であることを考えると、式(4)は刻み幅 1 で関数を離散化していると考えることができる。これを、任意の刻み幅 Δ𝑥 で離散化するため、格子点 𝑗(∈ ℤ) における、𝑓(𝑥) の値を、改めて 𝑓(𝑗Δ𝑥) ≡ 𝑓𝑗とおき、近似係数 𝑠0,𝑘 と等しいとみなすと、
𝑥𝑓(𝑥) = ∑ 𝑓𝑗 𝜙0,𝑗 ( )
Δ𝑥

(6)

𝑗∈ℤ

となる。
次に、Wavelet Galerkin 法による離散化式を、1 次元線形移流方程式を例に導出する。未知のスカラー関数 𝑢(𝑥, 𝑡) に関する 1 次元線形移流方程式は以下のようにあらわされる。
𝜕𝑢
𝜕𝑢
(7)
+𝑐
=0
𝜕𝑡
𝜕𝑥
ここで、𝑐 は移流速度を表す。この方程式の一般解は、(8)𝑢(𝑥, 𝑡) = 𝑢(𝑥 − 𝑐𝑡)で表される任意の関数であり、これは、ある時刻での分布が、𝑡 秒後に 𝑐𝑡 だけ平行移動する
ことを表す。重み付き残差法は、(7)式において、未知の関数 𝑢 を、𝑢̂ と近似した時に生じる
残差
𝜕𝑢̂
𝜕𝑢̂
(9)
+𝑐
=𝑟
𝜕𝑡
𝜕𝑥
に、任意の重み関数 𝑤 をかけて積分したものがゼロとなる、つまり
∞
𝜕𝑢̂
𝜕𝑢̂
∫ 𝑤 ( + 𝑐 ) d𝑥 = 0
(10)
𝜕𝑡
𝜕𝑥
−∞
となるように 𝑢̂ を定めるという方法である。さらに、(10)式において、近似関数と重み関数を同じ関数形で表す方法を Galerkin 法という。その中でも重み関数に区分多項式を用いて積分を各要素で実行し、格子点での関数値に関する連立一次方程式を導く方法を有限要素法といい、重み関数に直交関数を用いることで関数を周波数の異なるそれぞれのモードに分解し、各基底関数の係数を求める方法をスペクトル法という。Wavelet Galerkin 法は、(10)式に、(6)式のように Wavelet 基底を使用し、基底関数の直交性を利用して、各格子点での関数値に関する常微分方程式を導く方法である。このように、Wavelet Galerkin 法は、基底関数の直交性が基本にあるため、スペクトル法に近い考え方の方法であると考えることができる。ここで、(6)式と同様に 𝑢 の近似関数を(10)式に代入すると、

4

d𝑢𝑗
d𝜙0,𝑗 𝑥𝑥𝑥
∫ 𝜙0,𝑘 ( ) ∑ {𝜙0,𝑗 ( ) + 𝑐𝑢𝑗( )} d𝑥 = 0
Δ𝑥
d𝑡
Δ𝑥
d𝑥 Δ𝑥

（11）

𝑗

⟺∑
𝑗

d𝑢𝑗
𝑥
𝑥
𝑥 d𝜙0,𝑗 𝑥∫ 𝜙0,𝑘 ( ) 𝜙0,𝑗 ( ) d𝑥 + ∑ 𝑐𝑢𝑗 ∫ 𝜙0,𝑘 ( )( ) d𝑥 = 0d𝑡Δ𝑥Δ𝑥Δ𝑥 d𝑥 Δ𝑥(12)𝑗

ここで、𝜈 = 𝑥/𝛥𝑥とおくと、∑𝑗d𝑢𝑗
1 d𝜙0,𝑗(𝜈)Δ𝑥d𝜈 = 0∫ 𝜙0,𝑘 (𝜈)𝜙0,𝑗 (𝜈)Δ𝑥d𝜈 + ∑ 𝑐𝑢𝑗 ∫ 𝜙0,𝑘 (𝜈)
d𝑡
Δ𝑥 d𝜈

(13)

𝑗

条件 1 の、スケーリング関数の正規直交性を考慮して左辺を整理すると 𝑘 = 𝑗 のときだけ、左辺第 1 項は値を持つから、
d𝜙0,𝑗
d𝑢𝑘
1
(𝜈)d𝜈 = 0+
∑ 𝑐𝑢𝑗 ∫ 𝜙0,𝑘 (𝜈)
d𝑡
Δ𝑥
d𝜈

(14)

𝑗

左辺第 2 項について、以下のように考える。
𝜙0,𝑘 (𝜈)d𝜙0,𝑗d𝜙(𝜈) = 𝜙(𝜈 − 𝑘)(𝜈 − 𝑗)
d𝜈
d𝜈

(15)

あらためて 𝑥 = 𝜈 − 𝑘 とおくと、∫ 𝜙(𝜈 − 𝑘)

d𝜙
d𝜙
d𝜙0,𝑙
(𝑥)d𝑥
(𝜈 − 𝑗)𝑑𝜈 = ∫ 𝜙(𝑥)(𝑥 − 𝑙)d𝑥 = ∫ 𝜙0,0 (𝑥)
d𝑥
d𝑥
d𝑥

(16)

ここで、𝑙 = 𝑗 − 𝑘 である。そこで、(14)式の左辺第 2 項を、(15)式、(16)式を用いて次のように変形する。
∑ 𝑐𝑢𝑗 ∫ 𝜙0,𝑘 (𝜈)
𝑗

d𝜙0,𝑗
d𝜙0,𝑙
(𝜈)d𝜈 = ∑ 𝑐𝑢𝑘+𝑙 ∫ 𝜙0,0 (𝑥)
(𝑥)d𝑥
d𝜈
d𝑥

(17)

𝑙

よって(14)式は以下のようになる。
d𝑢𝑘
1
d𝜙0,𝑙
(𝑥)d𝑥 = 0+
∑ 𝑐𝑢𝑘+𝑙 ∫ 𝜙0,0 (𝑥)
d𝑡
Δ𝑥
d𝑥

(18)

𝑙

さらに、スケーリング関数 𝜙(𝑥) の台の幅が 𝐿 であることから、𝑙 ≤ −𝐿 および 𝐿 ≤ 𝑙 のとき
∫ 𝜙0,0 (𝑥)d𝜙0,𝑙(𝑥)d𝑥 = 0d𝑥(19)となることを考えると、あらためて格子点 𝑗 に関して
𝐿−1

d𝑢𝑗
1
(1)
+
∑ 𝑐𝑢𝑗+𝑙 𝑎𝑙 = 0
d𝑡 Δ𝑥

(20)

d𝑞 𝜙0,𝑙
(𝑥)d𝑥
d𝑥 𝑞

(21)

𝑙=−𝐿+1を得る。ここで、(𝑞)𝑎𝑙≡ ∫ 𝜙0,0 (𝑥)5である。また、時間微分項を前進差分によって離散化すると、
𝑐Δ𝑡
𝑢𝑗𝑛+1 = 𝑢𝑗𝑛 −
Δ𝑥

𝐿−1
(1)

∑ 𝑢𝑗+𝑙 𝑎𝑙(22)𝑙=−𝐿+1となる。𝑞 はスケーリング関数の導関数の階数を表し、拡散方程式など、高次の導関数の項がある場合でも同様の手順で計算できる。
また、スケーリング関数の多くは、明示的な関数形が表されないため、今回の検討では、
数値計算ソフトのひとつである Matlab の、Wavelet Toolbox 内に組み込まれているデータを出力し、台形公式により数値積分を行った。式(22)のスキームを基本として、本研究のスキームを構築していく。例として、coiflet のスケーリング関数を図 2-1 に、式(21)の計算結果を図 22 に示す。ここで coif の後の数字は、コイフレット次数といい、これが大きいほど台が広く、関数は滑らかになる。

(1)

図 2-1 coiflet スケーリング関数 : 𝜙

図 2-2 基底関数 coiflet による展開係数 : 𝑎𝑙式(20)は、空間微分を、着目している格子点の前後の値に係数をかけて足すことで勾配を求めているとみることができ、これは差分法の離散化式に類似していると考えられる。例えば、2次精度の中心差分では、空間方向の勾配は
d𝑢𝑗
1 1
1
+𝑐
( 𝑢𝑗+1 + 0 × 𝑢𝑗 − 𝑢𝑗−1 ) = 0
d𝑥
Δ𝑥 2
2

(23)

と表すことができる。このときの展開係数も併せて図 2-2 に示した。図 2-2 のグラフが反対称な形をしている点を考えると、本手法による離散化は、空間方向に中心差分を用いるのと似ていると考えることができる。次節では、式(22)による移流方程式の計算結果を示す。

2.3

線形移流方程式による検証と課題点

前節で導出した離散化式を用いて、1 次元線形移流方程式の数値シミュレーションを行う。
表 2-1 に計算条件を示す。
基底関数は coif1 を用いた。初期条件としてサインカーブを入力し、境界条件として周期境界条件を課した。

6

表 2-1 移流方程式計算条件

0  x 1

計算領域
空間刻み幅 Δ𝑥
時間刻み幅 Δ𝑡
移流速度 𝑐0.005
0.00011

このとき、Courant 数 ν を
𝑐Δ𝑡
(24)
Δ𝑥
と定義すると、表 2-1 の条件では ν = 0.02 である。この問題は、計算領域の幅が 1、移流速度が 1、周期境界であるため、解は周期 𝑇 = 1 で同じ波形となるはずである。そこで、9 周期分計算した結果を図 2-3 に示す。このとき、計算が正しく行われていれば、𝑡 が整数のときはすべて初期条件と同じ波形になるはずである。
ν≡

図 2-3 Wavelet Galerkin 法によるサインカーブの移流計算(coif1)
図 2-3 を見ると、1 周期後の解の分布はほぼ初期分布と同じになっており、精度よく計算できているようにみえる。しかし、計算を進めると、数値振動が発生することが確認できる。

計算では、このあと振動が増幅し、やがて発散した。次に、初期分布として矩形波とよばれる、不連続点をもつ波形を与えた時の結果を図 2-4 に示す。

図 2-4 Wavelet Galerkin 法による矩形波の移流計算

7図 2-4 をみると、数値振動により、初期の波形を保つことができていないことがわかる。すな

わち前節で推測した通り、このスキームは差分法における空間中心差分と同様に発散することが確認された。この原因は、スキームの特性が、解くべき方程式の性質と合致していないためであると考えられる。流れの問題では、情報は上流から下流へ一方的に伝播するため、ある点の時間発展を計算するのに必要なのは、着目点からみて上流側の情報だけである。しかし、中心差分などのスキームでは、着目点の上流側と下流側の情報が等しく含まれているため、方程式の物理的な性質とスキームの特性が矛盾する。これにより、中心差分で移流方
程式を解こうとすると、計算が不安定になるのである。Wavelet Galerkin 法の計算結果も、同様の理由で発散したと考えられる。このような不安定性を除去するために、差分法や有限要素法では様々な方法が考案されている。上流の点のみを使って移流項を差分近似する風上差
分法や、時間微分項の離散化を高次化することなどが、その代表例である。次章では、それらの方法の一つである、Taylor Galerkin 法の考え方を導入する。

8

第3章

Wavelet Taylor Galerkin 法の概要とスキ

ームの検証
3.1

Taylor Galerkin 法の概要

前章で課題となった、スキームの安定性について、Taylor Galerkin 法の考え方により、改善を試みる。ある時刻の解 𝑢𝑗𝑛 = 𝑢(𝑗Δ𝑥, 𝑛Δ𝑡) が既知の時、次時刻ステップにおける解 𝑢𝑗𝑛+1 =𝑢(𝑗Δ𝑥, (𝑛 + 1)Δ𝑡) は、Taylor 展開により
𝜕𝑢 Δ𝑡 2 𝜕 2 𝑢++ 𝑂(𝛥𝑡 3 )𝜕𝑡
2 𝜕𝑡 2
ここで、式(7)で表される基礎方程式と、その 1 階時間微分より
𝜕𝑢
𝜕𝑢
= −𝑐
𝜕𝑡
𝜕𝑥
2
𝜕 𝑢
𝜕
𝜕𝑢
𝜕 𝜕𝑢
𝜕2𝑢
2
=
(−𝑐
)
=
−𝑐
=
𝑐
𝜕𝑡 2 𝜕𝑡
𝜕𝑥
𝜕𝑥 𝜕𝑡
𝜕𝑥 2
となるので、これらを式(25)に代入し、高次の微小項を無視すると、
𝑢𝑗𝑛+1 = 𝑢𝑗𝑛 + Δ𝑡

(25)

(26)
(27)

𝜕𝑢
Δ𝑡 2 𝜕 2 𝑢
(28)
+ 𝑐2
𝜕𝑥
2 𝜕𝑥 2となる。右辺第 3 項は拡散項に相当し、これが数値振動を抑制する効果を持つ。式(28)を、Galerkin 法で空間方向に離散化する方法を Taylor Galerkin 法といい、差分法で離散化する方法
を Lax Wendroff 法[11]という。ここでは、Wavelet Galerkin 法で離散化するので、この手法をWavelet Taylor Galerkin 法とよぶ。式(28)を Wavelet Galerkin 法で離散化すると
𝑢𝑗𝑛+1 = 𝑢𝑗𝑛 − 𝑐Δ𝑡𝑐Δ𝑡
𝑢𝑗𝑛+1 = 𝑢𝑗𝑛 −
Δ𝑥

𝐿−1

∑

𝐿−1

(1)
𝑛
𝑢𝑗+𝑙
𝑎𝑙 + 𝑐 2

𝑙=−𝐿+1Δ𝑡 2(2)𝑛∑ 𝑢𝑗+𝑙
𝑎𝑙
2Δ𝑥 2

(29)

𝑙=−𝐿+1となる。また、式(28)を 2 次精度の中心差分法で離散化した Lax-Wendroff 法は、
𝑐Δ𝑡 𝑛
Δ𝑡 2
𝑛
𝑛
(30)
(𝑢𝑗+1 − 𝑢𝑗−1
) + 𝑐2(𝑢𝑛 − 2𝑢𝑗𝑛 + 𝑢𝑗−1)
2Δ𝑥
2Δ𝑥 2 𝑗+1となる。次節では、式(29)を用いて矩形波の移流計算を行うことで、スキーム改良の効果を検証する。(29)式右辺の、2 階の展開係数は図 3-1 のようになる。
𝑢𝑗𝑛+1 = 𝑢𝑗𝑛 −図 3-1 基底関数 coiflet による 2 階の展開係数

9

3.2

Wavelet Taylor Galerkin 法による移流方程式の計算結果

前節で導出した Wavelet Taylor Galerkin 法のスキームを用いて、表 2-1 の条件で矩形波の移流計算を行う。coif1 から coif5 の基底関数を用いた場合の結果と、差分法の Lax Wendroff 法による結果を図 3-2 に示す。

(a) coif1(b) coif2(c) coif3(d) coif4(e) coif5(f) 差分法(Lax Wendroff 法)
図 3-2 Wavelet Taylor Galerkin 法による矩形波の移流計算
図 3-2 より、Wavelet Taylor Galerkin 法を用いることで、数値振動を抑制し、安定に計算でき

10

ていることが確認できる。また、同条件では、Lax Wendroff 法よりも初期の波形を精度よく保っている。ここで、それぞれの基底関数を用いた場合の、計算終了時の各格子点における厳密解との誤差の二乗和を求めたものを図 3-3 に示す。空間の分割数が同じであるので、誤差の二乗和が小さいものほど精度がよいといえる。

図 3-3 基底関数の種類と格子点における誤差の二乗和

今回の条件では、coif5 が最も精度が良いことが分かった。次に、本手法が時間発展につれて誤差が増幅しないか、また、増幅するとしたらそれはどのような条件のときであるかを調べる。スキームに対するこのような分析を、安定性解析といい、移流方程式のような線形のス
カラー方程式のスキームに対しては、Von Neuman の安定性解析[12]がよく用いられる。

Wavelet Taylor Galerkin 法に対する安定性解析

3.3

von Neumann の安定性解析は、線形でスカラーの方程式に対する差分スキームの安定性解
析における代表的な方法である。ここでいう安定性とは、偏微分方程式の解に含まれるどんな成分の波形も、スキームの時間発展にともなって振幅が増幅しないという性質をいう。まず、対象となる偏微分方程式のある時間ステップでの解が𝑢𝑗𝑛 = 𝑔𝑛 𝑒 𝑖𝑗𝜃(31)のようにあらわせるとする。本節では、𝑖 を虚数単位とする。式(31)は、偏微分方程式の Fourier
級数解の一成分と解釈することもできる。スキームが安定であるとは、式(31)を時間発展させたときに、どんな 𝜃 に対しても解の振幅が増幅しないことと同じである。式(31)を式(29)の右辺に代入し、Euler の公式などをつかって整理すると
𝐿−1
𝐿−1
1 2 (2)
(1)
𝑛+1
𝑛 𝑖𝑗𝜃
𝑢𝑗
= 𝑔 𝑒 {[1 + (−𝜈𝑎0 + 𝜈 𝑎0 ) + ∑ 𝐴𝑘 cos𝑘𝜃] − 𝑖 ∑ 𝐵𝑘 sin𝑘𝜃}
2
𝑘=1
𝑘=1

(32)

ただし
1
1
(1)
(2)
(1)
(2)
(33)
𝐴𝑘 ≡ (−𝜈𝑎−𝑘 + 𝜈 2 𝑎−𝑘 ) + (−𝜈𝑎𝑘 + 𝜈 2 𝑎𝑘 )
2
2
1
1
(1)
(2)
(1)
(2)
(34)
𝐵𝑘 ≡ (−𝜈𝑎−𝑘 + 𝜈 2 𝑎−𝑘 ) − (−𝜈𝑎𝑘 + 𝜈 2 𝑎𝑘 )22である。von Neumann の安定性解析では、増幅率の絶対値がどんな波数成分に関しても、1 を超えることがないとき、そのスキームは安定であるという。また、そのようになるクーラン数の範囲を安定条件とよぶ。式(32)において増幅率は以下のようになる。
11

2

𝐿−1

𝐿−1

2 1/2

1
(1)
(2)
𝑃 ≡ {[1 + (−𝜈𝑎0 + 𝜈 2 𝑎0 ) + ∑ 𝐴𝑘 cos𝑘𝜃 ] + (∑ 𝐵𝑘 sin𝑘𝜃) }
2
𝑘=1

(35)

𝑘=1

それぞれの基底関数に対して、波数とクーラン数をパラメータとして、数値的に安定性解析を行った結果を図 3-4 に示す。

(a) coif1(b) coif2(c) coif3(d) coif4(e) coif5
図 3-4 Wavelet Taylor Galerkin 法 (coiflet 基底)の安定性解析
12ここで、𝜃 = 0 は空間に関して一定の成分、𝜃 = 𝜋 は格子点ひとつおきに符号を変える波であり、離散点で再現できる最も高周波な成分である。図 3-4 において、どの 𝜃 においても増幅率Pが 1 を超えていないとき、そのクーラン数で数値解は安定であるといえる。たとえば、
coif1 では、図 3-5 のように ν = 0.35 では安定だが、ν = 0.4 では不安定となる。

(a) ν = 0.35(b) ν = 0.4
図 3-5 安定性解析の検証(coif1)

以上より、von Neumann の安定性解析が Wavelet Taylor Galerkin 法にも有効であると考えられる。よって、この安定性解析の結果を基に、各基底関数を用いた場合のスキームの性質を考察する。まず、coif1 とそれ以外の基底関数では、増幅率の変化の傾向が全く異なった。coif1
では、0.35 < ν < 0.4 の区間に明確な安定条件が存在し、発散に向かってゆくのは 𝜃 = 𝜋 付
近の波数成分であると推測される。一方で、coif2 から coif5 では、あるクーラン数で急激にスキームが不安定になるのではなく、クーラン数が大きくなるにつれて徐々に 𝜃 = 2𝜋/3 付
近で増幅率が大きくなり、スキームを不安定にすることがわかる。しかし、coif2 から coif5 の場合、増幅率が急激に 1 を大きく超えることはなく、少しずつ誤差が増幅して発散に向かっていくことが推測される。以上より、以降の検討では、時間刻みを十分細かくし、基底関数には coif5 を用いることにする。

3.4

Godunov の定理とスキームの非線形化

前節で、Taylor Galerkin 法の考え方を導入することにより、矩形波の移流計算の精度を向上させることができた。しかしながら、数値振動を完全に除去することはできなかった。数値
振動の発生は、スキームの精度と密接にかかわっており、それを表現しているのが Godunovの定理[12]である。Godunov の定理によると、移流方程式に対する 2 次精度以上のどのような線形スキームも、解の単調性を維持することができない。つまり、高次精度の線形スキーム
では必ず数値振動が発生することを意味する。Wavelet Taylor Galerkin 法は、式(29)のように、各格子点の数値解の線形結合により次時刻の数値解を求めるため、線形スキームである。よって、Wavelet Taylor Galerkin 法では、数値振動の発生が保証されてしまう。そこで、精度を落とさずに数値振動を抑制させる方法が、スキームを非線形にすることである。古典的な方法は、解の分布によって変化する係数を持つ人工粘性項を付加することである。そこで本節
では、さらに数値振動を抑制するために、文献[12]を参考に、人工的な粘性項を陽的に付加す13

ることを検討する。Wavelet Taylor Galerkin 法により、計算された次時刻での関数値を𝑢̅𝑗 とすると、𝑢̅𝑗 には数値振動が含まれる。そこで、次の時間ステップに進む前に、
𝑢𝑗𝑛+1 = 𝑢̅𝑗 + 𝜀|𝑢̅𝑗+1 − 2𝑢̅𝑗 + 𝑢̅𝑗−1 ||𝑢̅𝑗+1 + 2𝑢̅𝑗 + 𝑢̅𝑗−1 |(𝑢̅𝑗+1 − 2𝑢̅𝑗 + 𝑢̅𝑗−1 )(36)のような処理を行う。ここで 𝜀 は経験的に決められる値である。式(36)の右辺第 2 項は、拡散項の形をしており、変動の大きいところでのみ、その係数部分が大きくなる。このように、数値振動をなまらせる効果を期待して陽的に付加される項を陽的人工粘性項とよぶ。coif5 を用いた場合について、陽的人工粘性項を加えた場合について検討した結果を図 3-6 に示す。

(a) 𝜀 = 1 × 10−4(b) 𝜀 = 1 × 10−2
図 3-6 陽的人工粘性の付加

図 3-5 より、陽的人工粘性の付加により、数値振動が抑制されていることが確認できる。また、係数 𝜀 が大きいと、人工粘性による数値拡散により、解がなまることが確認できる。

14

第4章

Wavelet Taylor Galerkin 法による Euler 方

程式の数値シミュレーション
4.1

提案手法による 1 次元 Euler 方程式の離散化

Euler 方程式は非粘性圧縮性流体の運動方程式である。また、Euler 方程式は質量、運動量、エネルギーの保存則を表す方程式を成分に持つ非線形偏微分方程式系でもある。CFD スキー
ムを開発するときは、Euler 方程式が用いられる。これは、実用上重要な高レイノルズ数流れ
では、物体付近を除いて粘性の効果が小さいこと、Euler 方程式は数学的な性質がよく知られていることなどが理由にある。[12]まず、1 次元非粘性流れに対する Euler 方程式の概要を示す。𝜌 を流体の密度、𝑢 を流速とすると、質量保存則の微分形は
𝜕𝜌 𝜕
(37)
(𝜌𝑢) = 0+𝜕𝑡 𝜕𝑥となる。式(37)は、検査面に流入する質量流量と流出する質量流量の差が、検査面内の密度の時間変化となることを表している。次に、𝑝 を流体の圧力とすると、検査面内の運動量は検査面の両端の圧力差によっても変化するから、運動量保存則は微分形で𝜕𝜕(38)(𝜌𝑢) +(𝜌𝑢2 + 𝑝) = 0𝜕𝑡𝜕𝑥のようにかける。次に、流体の単位体積当たりの全エネルギーを 𝑒 とすると、エネルギー保存則の微分形は𝜕𝑒 𝜕(39){(𝑒 + 𝑝)𝑢} = 0+𝜕𝑡 𝜕𝑥となる。ここで、全エネルギーとは、熱力学的内部エネルギーと運動エネルギーの和を表す。
式(39)は、検査面内のエネルギーの変化は、検査面両端でのエネルギー流束の差と、圧力がする仕事の差によっておこることを表す。さらに、理想気体を仮定する場合、状態方程式
1(40)𝑝 = (𝛾 − 1) (𝑒 − 𝜌𝑢2 )2を連立することで、未知数と式の数が一致し、解くことができる。式(37)、式(38)、式(39)をまとめて Euler 方程式とよび、以下のようにあらわす。
𝜕𝑸 𝜕𝑬
(41)
+
=0
𝜕𝑡 𝜕𝑥
ここで、𝑸 と 𝑬 は保存量、流束とよばれ、
𝜌𝜌
(42)
𝑸 = {𝜌𝑢} = {𝑚}
𝑒
𝑒
𝑚
𝜌𝑢
1
(𝛾 − 1)𝑒 + (𝛾 − 3)𝑚2 /𝜌2(43)
𝑬 = { 𝑝 + 𝜌𝑢2 } =1(𝑒 + 𝑝)𝑢32𝛾𝑒𝑚/𝜌 − (𝛾 − 1)𝑚 /𝜌
{}2である。ここで、𝑚 = 𝜌𝑢 である。
Euler 方程式の Taylor Galerkin 法による離散化式を導くために、保存量 𝑸 の 2 次の Taylor
展開を考える。

15

𝜕𝑸 1 2 𝜕 2 𝑸
(44)
+ Δ𝑡
𝜕𝑡 2
𝜕𝑡 2
ここで、右辺第 3 項の、保存量の時間 2 階微分をどのように計算するかが問題となるが、ヤコビ行列
𝑸𝑛+1 = 𝑸𝑛 + Δ𝑡

𝐀≡

𝜕𝑬
=
𝜕𝑸

0
1
(𝛾 − 3)𝑚2 /𝜌
2
2
3
3
(𝛾
[−𝛾𝑚𝑒/𝜌 + − 1)𝑚 /𝜌10−(𝛾 − 3)𝑚/𝜌𝛾−1

3𝛾𝑒/𝜌 − (𝛾 − 1)𝑚2 /𝜌2
2

𝛾𝑚/𝜌]

(45)

を導入すると、∂2 𝑸 𝜕 𝜕𝑸
𝜕
𝜕𝑬
𝜕
𝜕𝑬
𝜕
𝜕𝑸
𝜕
𝜕𝑬
(46)
= ( ) = (− ) =(− ) = − (𝐀 ) =
(𝐀 )
2
∂𝑡
𝜕𝑡 𝜕𝑡
𝜕𝑡
𝜕𝑥
𝜕𝑥
𝜕𝑡
𝜕𝑥
𝜕𝑡
𝜕𝑥
𝜕𝑥
となり、既知の量だけで保存量の 2 階時間微分を求めることができる。式(46)を式(44)に代入
すると、
𝜕𝑬𝑛 1 2 𝜕
𝜕𝑬𝑛
(47)
𝑸𝑛+1 = 𝑸𝑛 − Δ𝑡
+ Δ𝑡
(𝐀𝑛
)
𝜕𝑥
2
𝜕𝑥
𝜕𝑥
となる。式(47)を、Wavelet Galerkin 法により、離散化する。例として、第 1 成分について考えると、
3

𝑄1𝑛+1 = 𝑄1𝑛 − Δ𝑡𝑛
𝜕𝐸1𝑛 1 2 𝜕
𝑛 𝜕𝐸𝑖
+ Δ𝑡
(∑ A1𝑖
)
𝜕𝑥
2
𝜕𝑥
𝜕𝑥
𝑖=1

3

= 𝑄1𝑛 − Δ𝑡

2 𝑛
𝜕𝐸1𝑛 1 2
𝜕A𝑛1𝑖 𝜕𝐸𝑖𝑛𝑛 𝜕 𝐸𝑖+ Δ𝑡 ∑ (
+ A1𝑖
)
𝜕𝑥
2
𝜕𝑥 𝜕𝑥
𝜕𝑥

(48)

𝑖=1

式(48)を Wavelet Galerkin 法で離散化すると、以下のようになる。
𝐿−1

𝑄1𝑛+1
= 𝑄1𝑛𝑗 −𝑗Δ𝑡(1)∑ 𝐸1𝑛 𝑗+𝑙 𝑎𝑙Δ𝑥𝑙=−𝐿+13

𝐿−1

1 Δ𝑡 2+∑( ∑2 Δ𝑥 2
𝑖=1

𝐿−1

∑

(1,1)

(0,2)

(A𝑛1𝑖 𝑗+𝑙 𝐸𝑖𝑛 𝑗+𝑚 𝑏𝑙𝑚 + An1𝑖 𝑗+𝑙 𝐸𝑖𝑛 𝑗+𝑚 𝑏𝑙𝑚 ))(49)

𝑙=−𝐿+1 𝑚=−𝐿+1

ただしd𝑞 𝜙0,𝑙d𝑟 𝜙0,𝑚
(50)
(𝑥)
(𝑥)d𝑥
d𝑥 𝑞
d𝑥 𝑟
である。また、式(49)中の 𝐸𝑖𝑛 𝑗 で、𝑛 は時間ステップ、𝑗 は格子点の番号、𝑖 は対応する保存則の方程式を表す。つまり、𝑖 = 1 ならば質量保存則における流束 𝜌𝑢、𝑖 = 2 ならば運動量保存則における流束 𝑝 + 𝜌𝑢2 、𝑖 = 3 ならばエネルギー保存則における流束 (𝑒 + 𝑝)𝑢 の、(𝑥, 𝑡) = (𝑗Δ𝑥, 𝑛Δ𝑡) における値を意味する。
(𝑞,𝑟)

𝑏𝑙𝑚

4.2

≡ ∫ 𝜙0,0 (𝑥)計算モデルの概要

本節では、今回の検討で用いた計算モデルについて示す。今回のモデルでは、1 次元衝撃波管問題とよばれるモデルを用いた。概念図を図 4-1 に示す。

16

図 4-1 1 次元衝撃波管

これは、管を二つに仕切り、一方に高圧高温の気体を、他方に低圧低温の気体を充てんし、瞬間的に仕切りを取り去ったときの流れを模擬している。仕切りが取り去られると、衝撃波が形成され、低圧気体中を伝播する。一方、高圧気体中では、有心膨張波が伝播し、高温気体と低温気体の接触面では、接触不連続面が形成される。粘性や熱伝導を考慮していないため、接触面で気体の混合や熱伝導が生じず、理論的には接触面も不連続なまま伝播する。表 4-1 のような Sod の条件[13]とよばれる初期条件を与えた時、𝑡 = 0.1 における理論解は図 4-2 のようになる。この問題を数値的に解く場合、接触不連続面が伝播とともになまらないか、衝撃波の近傍で数値振動がどの程度発生するか、といった点に注意する必要がある。
表 4-1 初期条件

高圧側初期圧力 𝑝𝐻

1.0

低圧側初期圧力 𝑝𝐿0.1

高圧側初期密度 𝜌𝐻1.0

低圧側初期密度 𝜌𝐿

0.125

比熱比 𝛾

1.4

図 4-2 1 次元衝撃波管問題の理論解 (𝑡 = 0.1)

4.3

計算結果

基底関数に coif5 を用いて以下の条件で計算を行った結果と、同条件での Lax Wendroff 法による結果を図 4-3 に示す。なお、両者とも式(36)の陽的人工粘性項を付加している。
17

表 4-2 計算条件

計算領域

0<𝑥<1

空間刻み幅 ∆𝑥0.001

時間刻み幅 ∆𝑡

0.00001

パラメータ 𝜀0.01

基底関数の種類

coif5(a) 計算領域全体(b) 衝撃波近傍
図 4-3 1 次元衝撃波管問題の計算結果

結果より、本手法による数値解は、Lax Wendroff 法による数値解に比べて数値振動は小さい
が、衝撃波の伝播速度が理論解よりも速く計算されることが確認された。

4.4

陽的人工粘性項の改良とその効果改善のために、本検討では、陽的に付加した人工粘性項の計算方法に着目した。簡単のため、スカラー方程式に戻り、式(36)を再掲する
𝑢𝑗𝑛+1 = 𝑢̅𝑗 + 𝜀|𝑢̅𝑗+1 − 2𝑢̅𝑗 + 𝑢̅𝑗−1 ||𝑢̅𝑗+1 + 2𝑢̅𝑗 + 𝑢̅𝑗−1 |(𝑢̅𝑗+1 − 2𝑢̅𝑗 + 𝑢̅𝑗−1 )(36)

ここで、右辺第 2 項が人工粘性項であり、これは 2 次精度の中心差分に似た式になっている。
そこで、Wavelet Taylor Galerkin 法においては、この人工粘性項も、Wavelet Galerkin 法における拡散項と同様の方法で計算されるべきであると考え、以下のように人工粘性の計算方法を
変更する。
(2)

𝑢𝑗𝑛+1 = 𝑢̅𝑗 + 𝜀

̅𝑗+𝑙 𝑎𝑙 ||∑𝐿−1
𝑙=−𝐿+1 𝑢

𝐿−1(2)∑ 𝑢̅𝑗+𝑙 𝑎𝑙(2)
𝐿−1|∑𝑙=−𝐿+1 𝑢̅𝑗+𝑙 |𝑎𝑙 || 𝑙=−𝐿+1(51)

式(51)の人工粘性を用いて計算を行った結果を図 4-4 に示す。比較のため、図 4-4(b)には、Lax
Wendroff 法による結果も示してある。

18

(a) 計算領域全体(b) 衝撃波近傍
図 4-4 人工粘性項の計算方法による比較結果より、衝撃波の伝播速度が、より理論解に近く、Lax Wendroff 法の結果と同等になっていることが確認できる。また、数値振動は明らかに Lax Wendroff 法よりも小さい。このことから、Wavelet Taylor Galerkin 法による Euler 方程式の数値シミュレーションで、数値解の安定性を確保するために人工粘性項を付加する場合、その計算には Wavelet Galerkin 法における拡散項と同じ計算方法を適用する必要があることがわかる。さらに、人工粘性の効果を決める経
験的なパラメータ 𝜀 について、次節で検討する。

4.5

人工粘性項の経験的なパラメータに関する検討本手法では、人工粘性項の計算において、経験的なパラメータ 𝜀 を導入した。本節では数値シミュレーションを複数行うことによって、𝜀 に対して数値解の精度がどのような傾向を示すのかを見出すことを試みた。図 4-5 に、前節と同条件で、𝜀 を変えて計算を行った時の、各格子点における数値解と理論解の誤差の二乗和を示す。縦軸の値が小さいほど数値解が理論解に近いことを表す。

図 4-5 各手法の精度の比較

19

図 4-5 より、Wavelet Taylor Galerkin 法を用いた場合では、改良した人工粘性項による数値解の方が、精度がよかった。また、Lax Wendroff 法との比較では、𝜀 = 0.05 程度の時、提案手法の方が精度がよかった。しかし、提案手法では、𝜀 < 0.05 で精度が悪化し、𝜀 が小さすぎると計算が発散した。そのため、詳細な安定条件の解明が今後の課題としてあげられる。また、
𝜀 の増減に対して、提案手法と Lax Wendroff 法では、精度の変化が異なる傾向を示したのも今後原因を追究する必要があると考えられる。

20

第5章
5.1

まとめ
得られた知見本研究では、任意の格子幅での離散化が可能な Wavelet Taylor Galerkin 法のアルゴリズムを開発し、線形移流方程式および Euler 方程式の数値シミュレーションを行った。その過程で、陽的に付加する人工粘性項に着目し、その計算方法を改良した。本研究で得られた知見を以下に示す。
1. 開発した Wavelet Taylor Galerkin 法のスキームは、クーラン数が十分小さい時は、矩形波のよ
うな不連続面を含む問題でも精度よく計算を行うことができる。
2. coiflet 基底を用いる場合は、時間刻みを十分小さくしたうえで、高次の基底関数を用いる方が精度がよい。
3. 開発したスキームで Euler 方程式の 1 次元衝撃波管問題のシミュレーションを行った結果、差分法と同じ陽的人工粘性項付加すると、衝撃波の伝播速度が理論解よりも大きく計算されることが分かった。
4. 陽的人工粘性項を Wavelet Galerkin 法の展開係数を用いて計算するスキームでは、衝撃波の伝播速度がより理論解に近くなった。
5. 人工粘性項の経験的なパラメータを変えて複数回シミュレーションを行い、精度を比較した
結果、差分法とは異なる傾向を示し、特定の条件では提案手法の方が精度がよいことを確認し
た。

5.2

今後の課題

本研究で明らかになった課題および、本研究で実施できなかった未検討の事項についてま
とめる。
1.
2.
3.
4.

パラメータ 𝜀 の合理的な設定方法を見出す。
計算時間や使用メモリなどの改善。
粘性項や熱伝導項の導入。
2 次元および 3 次元問題への拡張。

1 については付録 D において予備的な検討を行っている。また 2 については、付録 C に示す
ような、多段解法の導入が効果的ではないかと考える。これは、1 ステップを数段に分けて時間進行させる方法で、1 ステップあたりの計算時間は長くなるが、時間刻みを長くとれるというメリットがある。これらの項目を解決後、3 および 4 が解決されることで、より実用的な流
れ場に本手法が適用可能になると考える。

21

謝辞
本研究を進めるにあたり、課題設定から論文執筆に至るまで非常に丁寧なご指導をいただ
きました山本亨輔助教に深く感謝いたします。また、ゼミで貴重なご助言をいただきました
松島亘志教授、学会などの機会にご指導をいただいた金川哲也助教に心より感謝いたします。
最後に、大学院での研究生活を支えてくれた研究室の皆様ならびに家族に心から御礼を申し
上げます。

22

参考文献
[1] 稲場智, 肖鋒, “THINC 法による圧縮性流体における不連続面の高解像度スキーム,” 第 29
回数値流体力学シンポジウム講演概要集(2017), 講演番号 C04-4[2] 青木一真, 牛島省, 鳥生大祐, “圧縮性流体に関する保存形陰的計算アルゴリズムの提案,”
土木学会論文集 A2(応用力学), 69 (2) (2013), pp.51-58.
[3] Liandrat, J. and Tchamitchian, Ph., “Resolution of the 1D regularized Burgers equation using aspatial wavelet approximation,” ICASE Report, 90-83 (1990).
[4] Donea, J., “A Taylor-Galerkin method for convective transport problems,” International Journal for
Numerical Methods in Engineering, 20 (1984), pp.101-119.
[5] Kumar, B. R. and Mehra, M., “A Wavelet-Taylor Galerkin method for parabolic and hyperbolic
partial differential equations,” International Journal of Computational Methods, 2 (1) (2005), pp.75-97.
[6] Kumar, B. R. and Mehra, M., “Wavelet-Taylor Galerkin method for the Buegers Equation,” BIT
Numerical Mathmatics, 45 (2005), pp.543-560.
[7] 山口貴史, 近藤弘一, “ウェーブレット級数展開に基づく数値シミュレーション,” 応用力学研究所研究集会報告 非線形波動及び非線形力学系の現象と数理(2006) , Article No. 35.
[8] 浜井祐治, 大田英輔, “Wavelet 分解を用いた Navier-Stokes 方程式解法スキーム構築の試み,”
日本機械学会論文集 B 編, 69 (2003), pp.249-257.
[9] 石岡圭一, スペクトル法による数値計算入門 (東大出版, 2004).
[10] 日本応用数理学会(編), 応用のためのウェーブレット (共立出版, 2016).
[11] Peter D. Lax, Burton Wendroff, “Systems of Conservation Laws,” Communications on Pure andApplied Mathematics, 17 (1964), pp.381-398.
[12] 藤井孝蔵, 流体力学の数値計算法, (東大出版, 1994).
[13] Gary A. Sod, “A survey of Several Finite Difference Methods for Systems of Nonlinear Hyperbolic
Conservation laws,” Journal of Computational Physics, 27 (1978), pp.1-31.

23

付録 A Wavelet Taylor Galerkin 法が 2 次以上の精度を持つことの検証
3.4 節で、Godunov の定理に基づき、Wavelet Taylor Galerkin 法のスキームに非線形な人工粘性項を導入した。そのとき、Godunov の定理の適用条件であるスキームの線形性には言及したが、スキームが 2 次精度以上であることは厳密には示していない。時間方向の 2 次の Taylor
展開に基づいたスキームであるため、Wavelet Taylor Glaerkin 法は 2 次精度ではあると考えられるが、時間の近似精度と、空間の近似精度は一般には分けて考えることはできない。そこ
で本節で、Wavelet Taylor Galerkin 法が時間・空間全体として 2 次以上の精度を持つことを数値的に確かめる。
Wavelet Taylor Galerkin 法の離散化式
𝐿−1

𝑢𝑗𝑛+1 = 𝑢𝑗𝑛 − 𝜈

∑

𝜈
(1)
𝑛
𝑢𝑗+𝑙
𝑎𝑙 +

2

𝑙=−𝐿+1

𝐿−1

2

(2)

𝑛
∑ 𝑢𝑗+𝑙𝑎𝑙(29)𝑙=−𝐿+1に Taylor 展開を適用し、何次の項まで消えるかを数値的に確認する。式(29)の左辺を時間、右辺を空間で Taylor 展開する(ν = cΔt/Δx)𝜕𝑢
1 𝜕2𝑢 2𝑛𝑢𝑗 +Δ𝑡 +Δ𝑡 + 𝑂(Δ𝑡 3 )
2
𝜕𝑡

2 𝜕𝑡

𝐿−1

= 𝑢𝑗𝑛 +∑

𝜈(1)[(−𝜈𝑎𝑙 +22

𝑙=−𝐿+1

𝜕𝑢
1 𝜕2𝑢
(2)
(𝑙Δ𝑥)2 + 𝑂(Δ𝑥 3 ))]𝑎𝑙 ) (𝑢𝑗𝑛 +
𝑙Δ𝑥 +
2
𝜕𝑥

A-1

2 𝜕𝑥

まず、0 次の項がキャンセルされるとき
𝐿−1(1)

𝑢𝑗𝑛 = 𝑢𝑗𝑛 +∑ [(−𝜈𝑎𝑙

+𝑙=−𝐿+1𝜈 2 (2) 𝑛𝑎 ) 𝑢𝑗 ]
2 𝑙

A-2

より
𝐿−1

−ν𝑢𝑗𝑛

∑

𝜈
(1)
𝑎𝑙 +

2

2

𝑙=−𝐿+1

𝐿−1
(2)

𝑢𝑗𝑛

∑ 𝑎𝑙

=0

A-3

𝑙=−𝐿+1A-3 が満たされる十分条件は
𝐿−1

∑

𝐿−1
(1)
𝑎𝑙 =

𝑙=−𝐿+1

(2)

∑ 𝑎𝑙

=0

A-4

𝑙=−𝐿+1であるが、coif5 の場合に A-4 の項を数値計算すると、数値的にほぼゼロになったので、Wavelet
Taylor Galerkin 法は A-2 を満たす。
次に 1 次の項は、
𝐿−1

𝜕𝑢
𝜕𝑢
𝜈 2 (2) 𝜕𝑢(1)
Δ𝑡 = −𝑐Δ𝑡 = ∑ [(−𝜈𝑎𝑙 + 𝑎𝑙 ) 𝑙Δ𝑥]
𝜕𝑡
𝜕𝑥
2
𝜕𝑥
𝑙=−𝐿+1

を満たす必要があるが、これはスキームが基礎方程式
24

A-5

𝜕𝑢
𝜕𝑢
+𝑐
=0
𝜕𝑡
𝜕𝑥
に適合しているための必要条件でもある。A-5 より
𝐿−1

𝜈 (1 − ∑𝑙=−𝐿+1

𝜈2(1)𝑙𝑎𝑙 ) +
2

A-6

𝐿−1
(2)

∑ 𝑙𝑎𝑙

=0

A-7

𝑙=−𝐿+1A-7 を満たすための十分条件は
𝐿−1
(1)

=1

A-8

(2)

=0

A-9

∑ 𝑙𝑎𝑙
𝑙=−𝐿+1

かつ
𝐿−1

∑ 𝑙𝑎𝑙
𝑙=−𝐿+1だが、coif5 の場合について計算すると、A-8、A-9 をほぼ満たしていた。次に、2 次の項について考える。
𝐿−1

1 𝜕2𝑢 2𝜈 2 (2) 1 𝜕 2 𝑢(1)(𝑙Δ𝑥)2 ]
Δ𝑡
=
∑
[(−𝜈𝑎
+
𝑎 )
𝑙
2 𝜕𝑡 22 𝑙
2 𝜕𝑥 2A-10

𝑙=−𝐿+1

また、基礎式を変形して
𝜕2𝑢
𝜕2𝑢
2
=𝑐
𝜕𝑡 2
𝜕𝑥 2

A-11

だから A-10 に代入して
𝐿−1

1 2
𝜈 2 (2)(1)𝜈 = ∑ [(−𝜈𝑎𝑙 + 𝑎𝑙 ) 𝑙 2 ]
2
2

A-12

𝑙=−𝐿+1

つまり
𝐿−1

𝐿−1

𝑙=−𝐿+1𝑙=−𝐿+1

1 2
1
(2)
(1)
𝜈 [1 −∑ 𝑎𝑙 𝑙 2 ] + ν ∑ 𝑎𝑙 𝑙 2 = 0
2
2

A-12

A12 を満たすための十分条件は
𝐿−1
(1)

A-13

(2)

A-14

∑ 𝑎𝑙 𝑙 2 = 0
𝑙=−𝐿+1かつ𝐿−1∑ 𝑎𝑙 𝑙 2 = 2
𝑙=−𝐿+1だが、coif5 について数値的に検証したところ、ほぼ A-13、A-14 を満たしていた。以上より、
coif5 を用いた Wavelet Taylor Galerkin 法のスキームは、2 次以上の精度を持つと考えられるので、本文にて Godunov の定理を適用した。
25

付録 B

Lax Wavelet Galerkin 法

第 2 章で、移流方程式に対する Galerkin 法のスキームが不安定であることを示し、第 3 章
では Taylor Galerkin 法の考え方に基づいて、安定なスキームを構築した。本節では、Taylor
Galerkin 法とは異なる方法で安定なスキームを構築する。まず、移流方程式に対する FTCS ス
キームの離散化式は∆𝑡𝑛B-1
𝑢𝑗𝑛+1 = 𝑢𝑗𝑛 − 𝑐(𝑢𝑛 − 𝑢𝑗−1)
2∆𝑥 𝑗+1のように表される。ここで、右辺第 1 項を、格子点 𝑗 の両隣の値の平均で置き換える。つま
り
𝑛
𝑛
𝑢𝑗+1
− 𝑢𝑗−1∆𝑡𝑛
B-2
𝑢𝑗𝑛+1 =−𝑐(𝑢𝑛 − 𝑢𝑗−1)
2
2∆𝑥 𝑗+1のように表す。このようなスキームを Lax 法とよぶ。この考え方を Wavelet Galerkin 法のスキームに導入すると、第 2 章の式(21)𝑐Δ𝑡
𝑢𝑗𝑛+1 = 𝑢𝑗𝑛 −
Δ𝑥

𝐿−1
(1)

∑ 𝑢𝑗+𝑙 𝑎𝑙

=0

(21)

𝑙=−𝐿+1

は以下のようになる。ここでは、このようなスキームを Lax Wavelet Galerkin 法とよぶこととする。
𝑢𝑗𝑛+1 =𝑛𝑛
𝑢𝑗+1
− 𝑢𝑗−1

2

𝐿−1

𝑐Δ𝑡
(1)
−
∑ 𝑢𝑗+𝑙 𝑎𝑙 = 0
Δ𝑥B-3

𝑙=−𝐿+1

Lax Wavelet Galerkin 法による移流方程式の計算結果を図 A-1 に示す。基底関数には coif5 を用いた。

(a) ν = 0.02(b) ν = 0.2
図 B-1 Lax Wavelet Galerkin 法による移流方程式の数値シミュレーション

Lax 法は数値拡散の効果が強いため、クーラン数が小さいときは分布がほぼ平坦になってしまう。そのため適用条件が限定される。

26

付録 C 多段階 Wavelet Galerkin 法
第 3 章で、Taylor Galerkin 法の考えに基づき、移流方程式に対する安定なスキームを構築した。この方法は、高次の Taylor 級数を考慮する方法であり、言い換えれば、時間方向の近似
を高次精度にすることが、スキームを安定にすると考えることができる。別の方法で、時間方向の近似を高次精度にする方法として、Runge Kutta 法に代表される多段解法がある。これは、次時刻の値を求めるために複数の段階に分けて計算する方法で、段数が多いほど高次精度にすることができる。たとえば、時間に関する常微分方程式d𝑢
C-1= 𝑓(𝑢(𝑡))dtという式に対して、次時刻の値を求めるために 3 段 3 次の公式を用いると
C-2
𝑢(1) = 𝑢𝑛 + Δ𝑡𝑓(𝑢𝑛 )11C-3
𝑢(2) = 𝑢𝑛 + Δ𝑡𝑓(𝑢𝑛 ) + Δ𝑡𝑓(𝑢(1) )
4
4
1
1
2
C-4
𝑢𝑛+1 = 𝑢𝑛 + Δ𝑡𝑓(𝑢𝑛 ) + Δ𝑡𝑓(𝑢(1) ) + Δ𝑡𝑓(𝑢(2) )663のように時間発展する。第 2 章のこれを式(20)
𝐿−1

d𝑢𝑗
1
(1)
+
∑ 𝑐𝑢𝑗+𝑙 𝑎𝑙 = 0
d𝑡 Δ𝑥

(20)

𝑙=−𝐿+1と組み合わせることで、3 段階 Wavelet Galerkin 法のスキームが構築できる。このスキームによる矩形波の移流計算の結果を図 B-1 に示す。基底関数には coif5 を用いた。

(c) ν = 0.02(d) ν = 0.2
図 C-1 3 段階 Wavelet Galerkin 法による移流方程式の数値シミュレーション

図より、クーラン数を大きくとっても計算が安定していることがわかる。一方で、この方法
は、計算を 1 時間ステップ進めるのに通常のスキームの 3 倍の計算量を要する。そのため、
時間ステップを長くとれるメリットと、1 時間ステップにかかる計算量のデメリットのバラ
ンスに注意する必要がある。

27

付録 D

人工粘性項のパラメータに関する追加検討
4.5 節で、人工粘性項に含まれる経験的なパラメータ 𝜀 に関して検討を行い、その決定方
法を今後の課題として挙げた。ここでは、計算条件を変化させたときに、最適な 𝜀 の値がどのような変化の傾向を示すか見出すために、予備的な検討を行った結果を示す。流れの問題
では、計算の質を議論する場合、クーラン数に着目することが多い。クーラン数とは、波や流れの物理的な伝播速度と、差分解の伝播速度の比であり、１次元の線形移流方程式に対しては、第 2 章で導入したように以下の式で定義される。
𝑐Δ𝑡
(24)
ν=
Δ𝑥
数値解の精度や安定性に関する議論は、クーラン数に基づいて議論されることが多い。ここで 𝑐 は解くべき偏微分方程式に含まれる物理的な波の速度であり、離散化によって変化することはない。離散化する場合に、計算者が決定できるのは格子幅 Δ𝑥 と時間刻み Δ𝑡 である。
そのため、4.5 節と同じ初期条件に対して、 Δ𝑡/Δ𝑥 をパラメータとして計算を行い、誤差の𝜀 依存性を確認した。
まず、格子幅と時間刻みを共に 2 倍にして計算を行った場合（ケース 1 とする）と、格子
幅のみ 2 倍にして計算を行った場合（ケース 2 とする）の、𝜀 と誤差の関係を図 D-1 に示す。

図 D-1 誤差の 𝜀 依存性（Δ𝑡/Δ𝑥 = 0.01, 0.005）

ケース 1 は、4.5 節と同様に 𝜀 = 0.005 付近で誤差が最小になっている。一方で、ケース 2 はグラフが全体的に左、つまり 𝜀 が小さくなる方向にシフトしている。クーラン数は小さいほうが計算が安定になることが多いことを考慮するとこの結果は、Δ𝑡/Δ𝑥 が小さい、すなわち
時間刻みが格子幅に比べて小さい場合、𝜀 を小さくとらないと、過剰な人工粘性が付加され28

ることを意味すると考えられる。以上より、最適な 𝜀 が、Δ𝑡/Δ𝑥 の関数で表される可能性が示唆されると考えている。

29