GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

2020/2/1 GFD オンラインセミナー

Wigner変換を⽤いた流体波動の局所分離解析
話題提供者: ⼤貫陽平 (九州⼤学 応⽤⼒学研究所)

論⽂: Onuki 2020, Quasi-local method of wave decomposition in a slowly varying
medium, Journal of Fluid Mechanics, 883, A56
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

1

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

⽬次

1. 導⼊: 地球流体における線形波動論の基礎
2. Wigner変換と擬微分作⽤素3. Wigner変換を⽤いた流体⽅程式の理論解析
4. シミュレーションデータ分析への利⽤

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

2

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

1. 導⼊: 地球流体における線形波動論の基礎
Fourier変換に基づく分散関係, 偏光関係の解析
空間に依存した系におけるWKB近似

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

3

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

はじめに:

回転成層流体である, ⼤気や海洋には多種多様な波動が存在する.
例: ⾳波, 慣性重⼒波, Rossby波, Kelvin波.
⼤気/海洋に⾒られる流体現象は波動⼒学の観点で説明されることが多い.
代表的なものとしては, 波動-平均流相互作⽤による惑星規模循環の形成.
波動間相互作⽤(海洋における内部重⼒波, ⼤気Rossby波)も重要なテーマ.
相互作⽤を考える際には, ⾮線形性が本質的に重要.
⼤抵の場合は摂動法に基づいて議論される.
線形波動論がすべての基礎になる.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

4

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

GFDにおける線形波動論の基礎

回転系の浅⽔⽅程式:
流速を(u, v), ⽔⾯変位をη, 重⼒加速度をg, コリオリパラメータをf , 静⽌状態での⽔深をH として,∂t u − f v = −g∂x η∂t v + f u = −g∂y η∂t η + ∂x (Hu) + ∂y (Hv) = 0

境界条件: 無限遠で, (u, v, η) → 0
この⽅程式にはどのような性質の波動が含まれるか?
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

5

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

初等的な考察: (f , H) が定数の場合

指数関数解: (u, v, η) = (u~, v~, η~)ei(kx+ly−ωt) を⽅程式に代⼊, 微分作⽤素を波数と周波数に置き換える:~ − f v~ = −igk η~−iω u~ = −iglη~−iω v~ + f u~ + ilH v~ = 0−iω η~ + ikH u

固有値問題を解いて周波数が定まる: ω = 0, ± f 2 + gH(k2 + l2 ) ... 分散関係
それぞれの分散関係に対応して, 従属変数の⽐である(u~, v~, η~)が定まる ... 偏光関係

独⽴した運動形態である, 慣性重⼒波モードと定常渦モードが得られる.

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

6

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

なぜ指数関数解を仮定して良い?
根底にあるのはFourier変換:

1~(k, l), v~(k, l), η~(k, l))ei(kx+ly) dkdl
(u, v, η) =∫ (u2π

それぞれの波数成分は独⽴に振る舞う.
あらゆる状態が波数成分の重ね合わせで表現できる.
与えられた(u, v, η)に対し, Fourier変換 → 波数空間で成分分解 → Fourier逆変換という⼿順を踏むことで, 独⽴したモード(慣性重⼒波と定常渦)を分離することができる.
実空間データ変換 波数空間データ 固有値分解, 射影

Fourier

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

分離データ(波数空間)

逆変換 分離データ(実空間)

Fourier

7

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析(f , H)が定数でない(x, yの関数となる)場合

波⻑に⽐べて f , H の変化の空間スケールが⼗分⼤きいと仮定.
⻑い空間スケールを導⼊し, μ(x, y) → (x, y), μ ≪ 1と置き換える.
空間微分が(∂x , ∂y ) → μ(∂x , ∂y ), と変更される.
解の形式を (u, v, η) = (u~(x, y), v~(x, y), η~(x, y))ei(θ(x,y)/μ−σt) と仮定して⽅程式に代⼊ (WKB近似), O(μ)の項まで考慮して, (u~, v~, η~, ∇θ, σ)の関係式を得る.
局所的な波数をk ≡ ∇θと定義すると, 分散関係式と偏光関係式が波数と空間位置に依存した形式で, σ = ω(x, k), (u~, v~, η~) ≡ u~ = u(x, k)と得られる.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

8

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

指数関数解とWKB近似解の⽐較
解の形式

分散関係式

偏光関係式

波動分離
指数関数解 u~ ei(k⋅x−σt) σ = ω(k) u~ = u(k) Fourier変換法
WKB近似解 u~ ei(θ(x)/μ−σt) σ = ω(x, k) u~ = u(x, k)?
(WKB近似における波数はk = ∇θ として定義)

Q1. WKB近似解に対応する波動分離⼿法はないか？
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

9

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

空間に依存した系に特有の性質

浅⽔⽅程式系において, 係数の空間変化(f (y), H(x, y))を考慮した場合:
惑星/地形性ベータ効果が⽣じて渦モードが有限の周波数を得る. (Rossby 波)
例: df /dy = β, H = const.とした時のRossby波:
−βk
2
ω=μ 2
+
O(μ
)
2
2
k + l + f /(gH)

浅⽔⽅程式で単に指数関数解を代⼊しただけでは導けない.
WKB近似を使って摂動項をゴリゴリ計算すれば導ける.
背景流の影響等を考慮するとかなり複雑な計算になる.

Q2. 分散関係式の摂動項を⾒通し良く求める⽅法はないか？

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

10

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

波動平均流相互作⽤研究の観点からエネルギー/擬運動量フラックスの解析. (群速度則: F = AC g を満たすように定式化)
e.g. Takaya & Nakamura (2001), Kinoshita et al. (2010), Kinoshita & Sato (2013)
基本場に対する擾乱について線形化された⽅程式の理論解析が主な内容.
WKB近似を前提とするが, 摂動項を無視することで実質的には指数関数解を仮定.
Kinosita & Sato (2013) は, プリミティブ⽅程式をvのみの式にしてから指数関数解を代⼊することで, Rossby波の分散関係式を導いている ... ややAdhocな印象.
問題設定をより複雑にして, 背景流の勾配とf の変化のスケールが同程度になるような場合(近慣性波やRossby波の分散関係が背景流シアを含む)には, 既存の⽅法論では太⼑打ちできない.

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

11

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

本研究: Wigner変換を⽤いて変数係数偏微分⽅程式の理論解析の⾒通しを良くする
複数種の波動を含む連⽴⽅程式を, 各波動成分が従う単独⽅程式に分解する.
連⽴偏微分⽅程式

Wigner変換

連⽴擬微分⽅程式

表象⾏列の
固有値分解

単独擬微分⽅程式

Fourier変換に基づく波動分離では代数⾏列の固有値分解であったのに対し, Wigner変換に基づく波動分離では表象⾏列の対⾓化を⾏う.

ポイント: 物理量データを固有波動成分へ射影する作⽤素を導く.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

12

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

2. Wigner変換と擬微分作⽤素
擬微分作⽤素と表象
スター積による計算
表象⾏列に関するアルゴリズム
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

13

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

通常のフーリエ変換では, 微分作⽤素を波数に置き換える.
∂x → ik

Wigner変換は, ⼀般の線形作⽤素を関数に変換する.

代表的な例: 量⼦⼒学のHamiltonianから古典⼒学のHamiltonianへの対応,
ℏ2 2
p2
−
∂x + V (x) →+ V (x)2m2m

上記の例では, ⼀座標xを残したまま, 微分作⽤素−iℏ∂x を運動量pに置き換えている.
Wigner変換によって, 変数の⾃由度は2倍になる: x → (x, p)
変換前の作⽤素と変換後の関数の関係は, 擬微分作⽤素とその表象の関係に対応.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

14

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

擬微分作⽤素(pseudo-diffential operator)は微分作⽤素を⼀般化した概念であり, フーリ
エ変換と逆変換を利⽤して, 積分変換の形式で
′
1
x
+
x
′ ip(x−x′ )/μ^F ϕ(x) =∫ f(, p) ϕ(x )edx′ dp2πμ2といった形で定義される. (上式でμは適当な定数であり, 量⼦論ではμ = ℏとする)
被積分関数に含まれるf (x, p)は, 擬微分作⽤素F^ の表象(symbol)と呼ばれる.
微分作⽤素も擬微分作⽤素の⼀種: f = (ip/μ)n ↔ F^ = ∂xn .
表象をf = (ip/μ)α , α ∈ Rという形に設定することで、逆ラプラシアンや分数
階微分なども擬微分作⽤素として定義できる.
Wigner変換とは, 擬微分作⽤素F^ をその表象f (x, p)へ対応させる変換である.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

15

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

Wigner変換の具体的⼿順:

対象となる作⽤素を積分変換の形式で表現する (Schwartzの核定理に基づく):F^ ϕ(x) = ∫ F (x, x′ )ϕ(x′ )dx′

核関数を相対位置座標についてFourier変換したものが作⽤素の表象になる:
x′x′−ipx′ /μf (x, p) = ∫ F (x + , x − ) edx′
22

例: 移流作⽤素 u(x)∂x を考えると, 核関数は F (x, x′ ) = u(x)δ′ (x − x′ )であり, これをフーリエ変換してf (x, p) = −i(p/μ)u + ∂x u/2を得る. 単に微分作⽤素を波数へ置き換えた場合に⽐べて追加の項が⽣じる.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

16

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

作⽤素の線形演算(⾜し引きやスカラー倍)は, 表象の演算に⾃然に対応する.
作⽤素の積は単純には対応付けられない.
作⽤素1 作⽤素2 作⽤素の線形演算
作⽤素の演算^^^作⽤素 F^αF^ + βGGF^ G
表象 f (x, p) g(x, p) αf (x, p) + βg(x, p) f (x, p) ⋆ g(x, p) (スター積)
積演算の定義が特殊なのは, 擬微分作⽤素が⼀般に⾮可換であるため基本となる２つの作⽤素をp^ ≡ −iμ∂x , x^ ≡ x と定義したとき, x^p^ = p^x^.
WKB近似ではμ ≪ 1を仮定. 最も粗い近似では, 作⽤素を可換とみなす. 作⽤素の⾮可換性によって摂動が⽣じる. このことがRossby波の発⽣に関与する. (f と∂y が⾮可換)
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

17

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

擬微分作⽤素の積に対応する表象の演算: スター積
擬微分作⽤素の定義式(p15)に⽴ち返って考察すると

f (x, p) ⋆ g(x, p)1′′′′′′ 2i(p′′ x′ −p′ x′′ )/μ
′
′′
′
′′
=
∬∬
f
(x
+
x
,
p
+
p
)g(x
+
x
,
p
+
p
)e
dx
dx
dp
dp
(πμ)2

が得られる. 被積分関数のf とgを(x, p)の周りでテイラー展開し, 積分を実⾏すると,
∞

(m+n)

(−1)n iμf (x, p) ⋆ g(x, p) = ∑
( )
m!n!
2
m=0,n=0∂ m+n f ∂ m+n g∂xm ∂pn ∂xn ∂pmとなる. 結果として, パラメータμのベキ級数展開で表現される. 以降は μ ≪ 1を仮定
し, ⾼次の項を摂動としてみなして話を進める.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

18

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

⾼次元系への拡張

⼀般の多変量関数 ϕ(x) : Rd → Cn に対して定義される擬微分作⽤素:
′
1
x
+
x
′ ip′ ⋅(x−x′ )/μ′^F ϕ(x) =
∫
f
(
,
p
)
ϕ(x
)e
dx
dp
d
(2πμ)
2

ここで, f はR2d で定義されたn × nの⾏列値関数. スター積は∣m∣+∣n∣(−1)∣n∣ iμf (x, p) ⋆ g(x, p) = ∑
( )
m!n!
2
m,n

n
n m
∇m
∇
f
(x,
p)∇
x
p
x ∇p g(x, p)(mとnには多重指数記法を⽤いている)以降は簡略化のため, 表象f (x, p)に対応する擬微分作⽤素をF^ = f (x^ , p^ )等と表記
する.

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

19

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

Wigner変換の視覚的表現

核関数の相対位置についてのFourier変換.
pは準局所的に定義された波数.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

20

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

Wigner変換のご利益

作⽤素の性質や微積演算を, 表象⾏列の性質や演算へ置きかえる.
エルミート作⽤素 ↔ エルミート⾏列
作⽤素の分解 ↔ 表象⾏列の分解
逆作⽤素の計算 ↔ スター積についての逆要素の計算
作⽤素の対⾓化 ↔ 表象⾏列の対⾓化

表象⾏列の対⾓化により, 流体⽅程式が含む波動成分を分離できる.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

21

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

逆作⽤素の計算

ある作⽤素 に対し, 以下の式を満たす逆作⽤素
^
L

−1
^
L

を求めたいとする:−1^^
LL = I( は単位⾏列). ここで と をWigner変換し, それぞれの表象をl(x, p), lI (x, p)とすると, 解くべき⽅程式は を既知として,
I^

−1
^
^
L L
l

となる. この式を解くため, 表象をl ⋆ lI = I

l = l0 + μl1 + μ2 l2 + …

と展開して代⼊する.

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

lI = lI0 + μlI1 + μ2 lI2 + …

22

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析(逆作⽤素続き) ⽅程式の各項をμのオーダーで分けると,
μ0 : l0 lI0 = I
1μ :l1 lI0 + l0 lI1 + l0 ⋆ lI0

1

=0

⋮

となる (表象aに対してa1 はμに⽐例した項を取り出すことを意味する). μ0 の式は通常の逆⾏列の計算により, lI0 = l−10 と解くことができ, ⾼次の式は
I
I
lI1 = −l−1
(
l
l
+
l
⋆
l
1 0
0
0
0

1

)

といった形で順次求めることができる. こうして得られたlI を表象にもつ擬微分作⽤素が所期の L−1 になる.

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

23

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

作⽤素の対⾓化与えられた作⽤素⾏列H = h(x^ , p^ )を, 対⾓形式に変換する問題を考える:−1^^^^ = diag(Ω^ 1, Ω
^ 2, ⋯ , Ω^ n)
U HU = Ω

表象⾏列に変換すると,⎛ ω1 (x, p)⎜0I⎜u⋆h⋆u =ω =⎜
⎜
⋮
⎝
0

0

…

ω2 (x, p)
…

⋱
0

⎞
⎟
⋮
⎟
⎟
⎟
0
ωn (x, p) ⎠0となる. 通常の⾏列の固有値問題に似ているが, 積がスター積に置き換わっている.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

24

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析(対⾓化続き1)
表象 h, ω, u, uI をμについて級数展開し, uとuI を以下のように書き換える:
1μ−1u = (I − u0 ⋆ u0 + μα) u0 + O(μ2 )
2
I

u

μ
−1
= u0 (I − u0 ⋆ u−1021− μα) + O(μ2 )解くべき⽅程式からμ0 とμ1 に⽐例した項を取り出すと, それぞれ
u0 h0 u−1
0 = ω0u0 ⋆ h ⋆ u−1
0

1

1
1
1
1
−1
−1
− u0 ⋆ u0 ω 0 − ω 0 u0 ⋆ u0 + αω 0 − ω 0 α = ω 122となる. 1つめの式は通常の⾏列の固有値問題であり、これは解けると仮定する.

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

25

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析(対⾓化続き2)μ1 の式の性質は, μ0 の⽅程式の結果によって異なる:(1) ω0 のすべての対⾓成分 (ω0,1 , ω0,2 , …) が異なる場合 (縮退がない)(2) ω0 の対⾓成分 (ω0,1 , ω0,2 , …) に重複がある場合 (縮退がある)(1) 縮退のない場合: ωとuのμ1 の補正項が次のように求まる:ω 1(k,k) = [u0 ⋆ h ⋆ u−1
0

1

]

(k,k)

− ω0,k [u0 ⋆ u−1
0

1

]

(k,k)

1
1
1
1 ω0,k + ω0,l−1
−1α(k,l) =[u0 ⋆ h ⋆ u0 ]−[u0 ⋆ u0 ]ω0,k − ω0,l2 ω0,k − ω0,l
(k,l)
(k,l)

(括弧書きの下付き添字は⾏列内における要素の位置を表す)

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

26

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析(対⾓化続き3)(2) 縮退がある場合:
μ1 の式は⼀般には計算できない.
→ 周波数が縮退した成分は分離することが出来ない.
特別な状況として縮退した固有値が0の場合 (ω0,i = ω0,i = … = 0):μ1 の式を計算することができ, 縮退が解ける.
→ 周波数0に縮退していた成分が空間の不均質性により, 固有の周波数を獲得する.
[具体的な計算⼿順については論⽂(Onuki 2020)を参照]
12

回転成層流体では, β 効果により定常渦モードがRossby波として周波数を獲得することに相当する.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

27

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

ここまでのまとめ:
Wigner変換は, ベクトル値関数に対する⼀般の作⽤素を, 表象⾏列に変換する.
作⽤素の積は表象⾏列のスター積に対応する.
逆作⽤素の計算や作⽤素の対⾓化等の演算は, 対応する表象⾏列を級数展開することで, 代数的な計算アルゴリズムに帰着される.

変数係数偏微分⽅程式の厄介な計算を, Wigner変換を通して⾒通しよく⾏うことができる.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

28

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

3. Wigner変換を⽤いた流体⽅程式の理論解析
流体⽅程式の⾮正準ハミルトン形式
線形化⽅程式の擬微分作⽤素表現単独⽅程式への分離
Wigner分布関数と輸送⽅程式
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

29

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

流体⽅程式の⼀般的表現

⾮粘性の流体運動を⽀配する⽅程式は, ⼀般に次の⾮正準ハミルトン形式で書かれる:δH^∂t υ = J (υ)
δυ

⋯ (⋆)

は状態ベクトルであり, 流速や密度等を表す.
H[v]はハミルトニアン汎関数で, エネルギーと適当なCasimirの和である.
†^^^ を満たす).
J はυ の関数として定義された歪エルミート作⽤素 (J = −J
地球流体⼒学における⾮正準ハミルトン形式の解説はShepherd (1990) あるいはSalmon (1998) を参照.
υ(x, t)

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

30

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

⽅程式の線形化とエネルギーの定義
がに代⼊,で極⼩となると仮定, そこからの擾乱をυ = Υ + υ′ と定義して⽅程式として⽅程式を線形化すると

H[υ] υ = Υ^ ≡ δ 2 H/δυ 2 , B^ = iA^J^ (Υ)A^(⋆)
A^ ∂t υ ′ = B^ υ ′ ⋯ (⋆⋆)
iAが得られる. ここでA^ は正値エルミート作⽤素, B^ はエルミート作⽤素.
この式から(擬)エネルギーの保存則が導かれる: E = (1/2) ∫ υ′† A^ υ′ dxとして,が成⽴.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

E ≥ 0; Ė = 0

31

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

作⽤素の分解前のページで定義した正値エルミート作⽤素A^ を, Wigner変換と表彰⾏列の演算(論⽂参照)によって, 以下のように分解する:†^^^
A = LL

ここで, は可逆作⽤素であり,ると, ⽅程式(⋆⋆)は
^
L

−1
^
L

が存在. 状態ベクトルを† ′^ψ=L υ/ 2と再定義す^ ψ ⋯ (⋆ ⋆ ⋆)i∂t ψ = Hと変形できる. ここでノルムとして,はエルミート作⽤素. 系のエネルギーは, ψのと書かれる.

−1
†−1
^
^
^
^
H = L BLE = ∫ ∣ψ∣2 dx

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

32

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

⽅程式の分離Wigner変換によってエルミート作⽤素に対応する表象⾏列 h(x, p) を定義, 固有値⽅程式u ⋆ h ⋆ uI = ωを解き対⾓⾏列 ω = diag(ω1 , ω2 , …) を求める. 状態変数を
ψ ′ = u(x^, p^ )ψ と再定義すると, ⽅程式(⋆ ⋆ ⋆)はi∂t ψ1′ = ω1 (x^, p^ )ψ1′i∂t ψ2′ = ω2 (x^, p^ )ψ2′⋮と分離できる. ここで, ωi (x^ , p^ )はそれぞれエルミート作⽤素であり, u(x^ , p^ )は, ユニ
タリ作⽤素である.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

33

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

分散関係と偏光関係前のページで得られた表象⾏列ω(x, p)の対⾓成分(ω1 , ω2 , …)はそれぞれの波
動成分の分散関係を表す.
対⾓化に⽤いた表象⾏列uの各列ベクトルが偏光関係を表す.
いずれの表象も, 背景パラメータの空間勾配(∇x f , ∇x H)と波束の⼤きさの有限性を考慮に⼊れて, 準局所的に定義された物理量と⾔える.
ユニタリ作⽤素u(x^ , p^ )は, 状態ベクトルψを固有空間へ写す射影作⽤素である.

以上の⼿順により, 空間に依存した系における波動分離を準局所的に⾏うことができる.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

34

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

密度作⽤素とWigner分布関数

ある波動成分の状態ベクトルψ(x)に対し, 次のように定まる作⽤素W^ を定義する:^ ϕ(x) ≡ ψ(x) ∫ ψ † (x′ )ϕ(x′ )dx′W

(ϕ(x)は任意の関数). このように定義されたW^ は, 量⼦⼒学では密度作⽤素と呼ばれる.
さらに, 密度作⽤素をWigner変換して得られる表象w(x, p)は, Wigner分布関数と呼ば
れる:
′
x′
x
†
−ip⋅x′ /μw(x, p) = ∫ ψ (x + ) ψ (x − ) edx′
22

Wigner分布関数は物理-波数空間内で定義されたエネルギー密度と解釈でき, 積分値は
全エネルギーに⼀致する: (1/2πμ)d ∬ w(x, p)dxdp = E .

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

35

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

輸送⽅程式

状態ベクトルの時間発展が i∂t ψ = ω(x^ , p^ )ψであるとき, Wigner分布関数の時間発展はi∂t w = ω ⋆ w − w ⋆ ωとなる. スター積をμで展開し, O(μ3 )の項を無視すると,∂t w + μ∇p ω ⋅ ∇x w − μ∇x ω ⋅ ∇p w = 0となる. この式は物理-波数空間内におけるエネルギー密度の伝達を記述するものであり, 輸送⽅程式(あるいは放射伝達⽅程式など)と呼ばれる.
輸送⽅程式の特性曲線は, ẋ = ∇p ω, ṗ = −∇x ω で定められる.
これは通常のRay-tracing ⽅程式である.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

36

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

エネルギーフラックス

輸送⽅程式を物理空間に射影(波数 p で積分)すると, エネルギー密度の時間発展がと得られる. ここで,

∂t E + ∇x ⋅ F = 0
1
E(x) =∫ w(x, p)dp
d
(2μπ)
μ
F (x) =∫ w(x, p)∇p ωdpd(2μπ)

である. エネルギーフラックス = エネルギー密度 × 群速度という形式になっている.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

37

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

ここまでのまとめ
⼀般の流体⽅程式に対し, Wigner変換を利⽤して波動成分を準局所的に分離する⼿順を構築.
Wigner分布関数を導⼊することで, 物理-波数空間内でのエネルギー伝達を記述.
群速度に⽐例した形式のフラックスを計算できることを⽰した.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

38

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

4. シミュレーションデータ分析への利⽤
浅⽔⽅程式の擬微分作⽤素表現数値モデルデータから波動成分を分離してエネルギーフラックスを描画

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

39

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

浅⽔⽅程式の擬微分作⽤素表現

流速を(u, v), ⽔深をh, 重⼒加速度をg, コリオリパラメータをf (y), ⽔底の地形をD(x, y)とする. ハミルトニアン(全エネルギー)を
hu2 + hv 2
g(h + D)2
H=∫ (
+
) dx
2
2

と定義して, 浅⽔⽅程式は⾮正準ハミルトン形式で次のように書ける:⎛ ∂t u ⎞ ⎛ 0
∂t v
=
−f /h
⎝ ∂t h ⎠ ⎝ −∂x

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

f /h
0
−∂y

−∂x ⎞ ⎛ δH/δu ⎞
−∂y
δH/δv
0 ⎠ ⎝ δH/δh ⎠40

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析(続き) 定常解として静⽌状態
u = 0, v = 0, h = H(x) ≡ H0 − D(x)を考え, 前述の⼿順(p31)に沿って擾乱について線形化された⽅程式を導く:
⎛ H
i
0
⎝ 0

0
H
0

0 ⎞ ⎛ ∂t u′ ⎞ ⎛00∂t v ′=−if H
g ⎠ ⎝ ∂t η ′ ⎠ ⎝ −ig∂x H

Hermite operator

状態ベクトルをψ = ( H/2u′ ,if H0−ig∂y H

−igH∂x ⎞ ⎛ u′ ⎞−igH∂yv′⎠ ⎝ η′ ⎠0Hermite operator

H/2v ′ ,g/2η ′ )と定義すると, 上の式はi∂t ψ = h(x
^, p
^ )ψ

と書ける. (表象⾏列hの具体的な形は次のページ)
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

41

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析(続き) 浅⽔⽅程式を特徴づけるエルミート作⽤素の表象⾏列h(x, p):
⎛

0
h=
−if
⎝ −cpx

if
0
−cpy
h0

⎛ 0
⎜
−cpx ⎞
⎜
−cpy
+ μ⎜
⎜ 0
⎜
0 ⎠
⎜ ic
x
⎝ −
2

ここで, p = (px , py ), c = gH , ∇x c = (cx , cy )とした.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

0
0
icy
−
2

icx ⎞
2 ⎟
icy ⎟
⎟
2 ⎟
⎟
⎟
0 ⎠

h1

42

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

浅⽔⽅程式の対⾓化定められた⼿順(p24-p27)に沿い, ⽅程式u ⋆ h(x, p) ⋆ uI = ωを解いて作⽤素を対⾓化する. 偏光関係を表す⾏列uは⎛ θpx − if py2θp⎜⎜ −θpx − if py
u=⎜
⎜
2θp
⎜
⎜
icpy
⎝
θ

θpy + if px
2θp−θpy + if px
2θp
−icpx
θ

となる. ( p = p2x + p2y , θ = f 2 + c2 p2 とした)
アルゴリズムの実⾏にはMaximaを利⽤.

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

cp ⎞
2θ ⎟
cp ⎟
⎟ + O(μ)
⎟
2θ ⎟
f ⎟
⎠
θ

43

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析(続き) 分散関係を表す ω = diag(ωG+ , ωG− , ωR )は,
2
2 2
βp
(2f
+
c
p ) f c(px cy − py cx )
x
2
2
2
2
ωG± = ± f + c p − μ [ 2 2
+
]
+
O(μ
)
2
2
2
2
2
2p (f + c p )f +c pc(−βcpx + 2f px cy − 2f py cx )2
ωR = μ
+
O(μ
)
2
2
2
f +c pとなる(β ≡ ∂y f とした). それぞれ慣性重⼒波とRossby波の周波数に対応する.
パラメータf , H の空間勾配を摂動として⽣じるRossby波の分散関係を正確に導くことができた.
β 効果により, 慣性重⼒波は波数ベクトルの⽅向によって分散関係が変化する.
(このことは最近 Perez et al. でも議論されている.)
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

44

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

数値モデルデータの分析

浅⽔系の数値シミュレーション.
⼆重周期境界条件の正⽅形領域で, ⽔底の形状を右上の図に設定.
初期の⽔⾯に変位を与えて, 波束の運動を観察する.
境界条件: ⼆重周期
格⼦点数: 128 × 128
格⼦幅∼変形半径/2⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

45

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

基本的な出⼒データ
左から: u, v, η

慣性重⼒波の速い伝搬と, 地形性Rossby波の遅い伝搬が⾒られる.

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

46

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

波動分離解析

射影作⽤素 u(x^ , p^ )を⽤いて波動成分を分離する.
左: 慣性重⼒波の信号 ψG
右: Rossby波の信号 ψR(上が実部, 下が虚部)(Rossby波の信号を解析する際, ωR が正の部分のみを取り出すフィルタ操作を⾏っている. 詳細は論⽂参照)
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

47

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

エネルギー解析

Wigner分布関数w(x, p)を計算, 物理空間に射影してエネルギー密度E(x)とフラックスF (x)を得る. [左:慣性重⼒波, 右:Rossby波]

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

48

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

エネルギー解析

ある瞬間のスナップショット.
それぞれの波の群速度に⽐例したフラックスを描画することができた.

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

49

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析

全体まとめ
流体波動を記述する偏微分⽅程式に対し, Wigner変換を⽤いることで, ⾒通
しよく⽅程式の変形操作を⾏うことができる.
Wigner変換により得られた表象⾏列の対⾓化により, 連⽴偏微分⽅程式に含まれる独⽴した波動成分が準局所的に分離される.
それぞれの波動の信号からWigner分布関数を定義し, 物理空間へ射影する
ことで, 群速度に⽐例したフラックスが得られる.
⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

50

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析(補⾜)

⼀般的な振動モードの分離
閉じた領域での固有値解析
離散モードとして得られる.
各モードは領域全体に広がった構造.
(右: Rossby-Haurwitz 波)
今回のテーマ無限に広い領域.
空間に局在した波束.

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

51

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析(補⾜) Maximaを⽤いた表象⾏列の対⾓化計算
define(th(x,y,px,py),sqrt(f^2+(c+d*y)^2*(px^2+py^2)))$
define(p(px,py),sqrt(px^2+py^2))$
define(h0(x,y,px,py),matrix([0,%i*f,(c+d*y)*px],[-%i*f,0,(c+d*y)*py],[(c+d*y)*px,(c+d*y)*py,0]))$define(h1(x,y,px,py),matrix([0,0,0],[0,0,%i*d/2],[0,-%i*d/2,0]))$
define(u0(x,y,px,py),matrix([(th(x,y,px,py)*px - %i*f*py)/(sqrt(2)*th(x,y,px,py)*p(px,py)),(th(x,y,px,py)*py + %i*f*px)/(sqrt(2)*th(x,y,px,py)*p(px,py)),(c+d*y)*p(px,py)/(sqrt(2)*th(x,y,px,py))],[(th(x,y,px,py)*px + %i*f*py)/(sqrt(2)*th(x,y,px,py)*p(px,py)),(th(x,y,px,py)*py - %i*f*px)/(sqrt(2)*th(x,y,px,py)*p(px,py)),-(c+d*y)*p(px,py)/(sqrt(2)*th(x,y,px,py))],[%i*(c+d*y)*py/th(x,y,px,py),-%i*(c+d*y)*px/th(x,y,px,py),f/th(x,y,px,py)]))$define(uc0(x,y,px,py),matrix([(th(x,y,px,py)*px + %i*f*py)/(sqrt(2)*th(x,y,px,py)*p(px,py)),(th(x,y,px,py)*px - %i*f*py)/(sqrt(2)*th(x,y,px,py)*p(px,py)),-%i*(c+d*y)*py/th(x,y,px,py)],[(th(x,y,px,py)*py - %i*f*px)/(sqrt(2)*th(x,y,px,py)*p(px,py)),(th(x,y,px,py)*py + %i*f*px)/(sqrt(2)*th(x,y,px,py)*p(px,py)),%i*(c+d*y)*px/th(x,y,px,py)],[ (c+d*y)*p(px,py)/(sqrt(2)*th(x,y,px,py)),-(c+d*y)*p(px,py)/(sqrt(2)*th(x,y,px,py)),f/th(x,y,px,py)]))$define(om(x,y,px,py),matrix([th(x,y,px,py), 0 , 0],[0, -th(x,y,px,py), 0],[0, 0, 0]))$

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp

52

GFD オンラインセミナー: Wigner変換を⽤いた流体波動の局所分離解析(続き)

define(du0dy(x,y,px,py), diff(u0(x,y,px,py) , y ) )$define(du0dpy(x,y,px,py) , diff(u0(x,y,px,py) , py) )$define(duc0dy(x,y,px,py) , diff(uc0(x,y,px,py) , y ) )$define(duc0dpy(x,y,px,py) , diff(uc0(x,y,px,py) , py) )$define(dh0dy(x,y,px,py), diff(h0(x,y,px,py), y ) )$
define(dh0dpy(x,y,px,py), diff(h0(x,y,px,py), py) )$
define(domdy(x,y,px,py), diff(om(x,y,px,py), y) )$
define(domdpy(x,y,px,py) , diff(om(x,y,px,py), py) )$
define(uhu1(x,px,py),u0(x,0,px,py).h1(x,0,px,py).uc0(x,0,px,py) +%i/2*(du0dy(x,0,px,py).dh0dpy(x,0,px,py).uc0(x,0,px,py)- du0dpy(x,0,px,py).dh0dy(x,0,px,py).uc0(x,0,px,py)+ u0(x,0,px,py).dh0dy(x,0,px,py).duc0dpy(x,0,px,py)- u0(x,0,px,py).dh0dpy(x,0,px,py).duc0dy(x,0,px,py)+ du0dy(x,0,px,py).h0(x,0,px,py).duc0dpy(x,0,px,py)- du0dpy(x,0,px,py).h0(x,0,px,py).duc0dy(x,0,px,py) ))$define(uu1(x,px,py),%i/2*(du0dy(x,0,px,py).duc0dpy(x,0,px,py)- du0dpy(x,0,px,py).duc0dy(x,0,px,py)))$define(a(x,px,py),uhu1(x,px,py)- om(x,0,px,py).uu1(x,px,py))$
factor(ratsimp(a(x,px,py)));define(lm(x,px,py),matrix([0,0,0],[0,0,0],[0,0,0]))$for k:1 thru 3 dofor l:1 thru 3 do
if k = l then lm(x,px,py)[k,l]:0
else lm(x,px,py)[k,l]:- 1/(om(x,0,px,py)[k,k] - om(x,0,px,py)[l,l]) * uhu1(x,px,py)[k,l]+ 1/2 * (om(x,0,px,py)[k,k] + om(x,0,px,py)[l,l]) / (om(x,0,px,py)[k,k] - om(x,0,px,py)[l,l])* uu1(x,px,py)[k,l]$define(u1(x,px,py),(- uu1(x,px,py)/2 + lm(x,px,py)).u0(x,0,px,py))$
factor(ratsimp(u1(x,px,py)));

⼤貫陽平 (九⼤応⼒研): onuki@riam.kyushu-u.ac.jp53