データ駆動型乱流研究に対する力学系的アプローチ

犬伏正信

東京理科大学 理学部応用数学科 / 大阪大学 大学院基礎工学研究科
GFDオンラインセミナー 第13回（2024年2月19日）

Van Dyke with Generative AI

2

Gen-2 by Runway

3

NG集

拡散⇨発散？

逆流？

異なる乱流？
4

阿波鳴門之風景（歌川広重）

Hiroshige Utagawa, The Whirlpools of Awa (1857) https://www.metmuseum.org/art/collection/search/45025

56データ駆動型乱流研究とその理論的基盤
AI: 原理不明のまま乱流の“シミュレーション”が可能
今後データから直接シミュレーションが可能に？
現状の多くのデータ駆動型研究は，モデル化性能の見極めや開発指針は経験やノウハウに依存．それらを保証する理論的基盤がない．

本研究：力学系理論に着目
力学系理論の概念や現象（Lyapunov次元やBlowout分岐）が最新のデータ駆動型乱流研究において重要になることを紹介．

7

おしながきカオスと乱流カオスとLyapunov指数，Lyapunov次元
Blowout分岐と横断Lyapunov指数
乱流の力学系的解釈

データ駆動型乱流研究
MLによる乱流モデルの低次元化限界（Lyapunov次元）[1]
乱流の小スケール渦の隷属性とデータ同化限界（横断Lyapunov指数）[2]まとめ[1] Linot and Graham, Dynamics of a data-driven low-dimensional model of turbulent minimal Couette ﬂow, J. Fluid Mech. 973, A42 (2023).
[2] Inubushi, Saiki, Kobayashi, and Goto, Characterizing small-scale dynamics of Navier-Stokes turbulence with transverse Lyapunov exponents: A dataassimilation approach, Phys. Rev. Lett. 131, 254001 (2023).

8カオスとLyapunov指数
カオスdx= f(x) (x(0) = x0 ∈ ℝN )の有界領域における決定論的力学系
dt
非周期軌道x(t) = φ t(x0)で初期値鋭敏性（λ1 > 0）のあるもの．

（最大）Lyapunov指数
誤差ベクトルδx(t)の指数的拡大率：∥δx(t)∥ ∝ ∥δx(0)∥e λ1t変分方程式∂f1∥δx(T )∥dδx= [Dfx]δx [Dfx]ij := i , (1 ≤ i, j ≤ N ) とするとき，λ1 = lim log．
dt
∂xj
T→∞ T
∥δx(0)∥()xN

φ t(x̃0)

x̃0x0

δx(0)

δx(t)

x2

φ t(x0)x1

9

Lyapunovスペクトルと次元
Lyapunovスペクトル{λ1, λ2, ⋯, λn} (λ1 ≥ λ2 ≥ ⋯ ≥ λn): λ1 + λ2 + ⋯ + λnはn次元体積の指数的増大率例）{λ1, λ2}: V2(t) = V2(0)exp((λ1 + λ2)t)，V2(t)は平行四辺形の面積

Lyapunov次元
Kaplan and Yorke (1979)：d = ñ +LxN
δx1(0)x0ñ∑j=1 λj| λñ+1 |ここで ñ := maxk|

k

λ ≥0
∑ j
j=1

λ1 + ⋯ + λñ ≥ 0, λ1 + ⋯ + λñ + λñ+1 < 0 （ ∴ Vñ+1(t) → 0 ）V0(t)δx2(0)

V2(t)Lorenz系

δx2(t)

δx1(t)x2z

φ t(x0)x1ydL ≃ 2.06 ( < 3)

xDoering and Gibbon (1995)

10Blowout分岐と横断Lyapunov指数
yExample [1-3]ℳ

平面内の1粒子系（x ∈ ℝ2）

xd 2x= − νvx − ∇V(x) + ( f0 sin ωt)ex2dt
ポテンシャル：V(x, y) = (1 − x 2)2 + y 2(x−p) + ky 4
(ν = 0.05, f0 = 2.3, ω = 3.5, k = 0.0025)
x· = v , v· = − νv + 4x(1 − x 2) + y 2 + f sin ωt{ y· = vy, v·y = − νvy − 2y(x − p) − 4ky 3xxx0V(x, y)p = − 1.9不変多様体ℳ := {(x, vx, y, vy, t) | y = vy = 0}
ℳ上のダイナミクスは“Duﬃng chaos”
x· = v , v· = − νv + 4x(1 − x 2) + f sin ωtxxx（※ pによらない）x(t)0x0 ∉ ℳを初期値とする解のダイナミクスは？

t[1] Sommerer and Ott, Nature (1993), [2] Ott et al., Phys. Rev. Lett. (1993), [3] Ott and Sommerer, Phys. Lett. A (1994).

11

Blowout分岐と横断Lyapunov指数
Blowout分岐[3]We use the term "blowout" because the feature of interest

is a dramatic departure of the system state from the vicinity of
the invariant manifold; the orbit is "blown out."p = − 2.0：ℳ内の“Duﬃng chaos”がアトラクタ
p = − 1.7：ℳ内の“Duﬃng chaos”はアトラクタではなくなる（chaotic saddle）

p = − 1.7x(t), y(t)p = − 2.0t
[1] Sommerer and Ott, Nature (1993), [2] Ott et al., Phys. Rev. Lett. (1993), [3] Ott and Sommerer, Phys. Lett. A (1994).

13

Blowout分岐と横断Lyapunov指数y

Blowout分岐の特徴付け[3]横断Lyapunov指数ℳに横断方向の摂動(δy, δvy)に対する安定性を定量化··変分方程式： δy = δvy, δvy = − νδvy − 2δy(x − p)（cf. y· = v , v· = − νv − 2y(x − p) − 4ky 3）yy(ℳxy

横断Lyapunov指数 λ⊥1δ(T )λ⊥ := lim log
T→∞ Tδ(0)δ(t) =δy[δy(t)]2 + [δvy(t)]2)

pc = − 1.7887でλ⊥の符号が変化
p < pcであればλ⊥ < 0，ℳは漸近安定[1] Sommerer and Ott, Nature (1993), [2] Ott et al., Phys. Rev. Lett. (1993), [3] Ott and Sommerer, Phys. Lett. A (1994).

14

乱流の力学系的解釈
乱流 ⊂ NSEによって定まる力学系のカオス
乱流のLyapunov指数： λ ∝ τ −1 (τ = (ν/ϵ)1/2: Kolmogorov time)[1]
近年の数値的研究（DNS）[2-4]：λτはわずかにRe数に依存して増大“we ﬁnd that it increases with the Reynolds number, but surprisingly faster than what ispredicted on dimensional grounds”[2][2][1] Ruelle, Phys. Lett. A 72 (1979)[2] Boﬀetta and Musacchio, PRL 119 (2017)[3] Mohan, Fitzsimmons, and Moser, PRF 2 (2017)[4] Arjun and Ho, PRL 120 (2018)15

（例）クエット乱流のLyapunov解析
平行平板間に起こる典型的な乱流（± U0で移動する上下平板によって駆動）
力学系の次元：d = 𝒪(105) (cf. Re = 400, MFU)

https://cns.gatech.edu/~gibson/research/index.html

16

（例）クエット乱流のLyapunov解析
平行平板間に起こる典型的な乱流（± U0で移動する上下平板によって駆動）
力学系の次元：d = 𝒪(105) (cf. Re = 400, MFU)
アトラクタ次元[1]：dL = 14.8 ≪ 𝒪(105)[1] M. Inubushi, S. Takehiro, and M. Yamada, Phys. Rev. E 92.2 (2015): 023022.

17おしながきカオスと乱流カオスとLyapunov指数，Lyapunov次元
Blowout分岐と横断Lyapunov指数
乱流の力学系的解釈

データ駆動型乱流研究
MLによる乱流モデルの低次元化限界（Lyapunov次元）[1]
乱流の小スケール渦の隷属性とデータ同化限界（横断Lyapunov指数）[2]まとめ[1] Linot and Graham, Dynamics of a data-driven low-dimensional model of turbulent minimal Couette ﬂow, J. Fluid Mech. 973, A42 (2023).
[2] Inubushi, Saiki, Kobayashi, and Goto, Characterizing small-scale dynamics of Navier-Stokes turbulence with transverse Lyapunov exponents: A dataassimilation approach, Phys. Rev. Lett. 131, 254001 (2023).

18

データ駆動型低次元モデル
ニューラルネット（autoencoder, Neural ODE, etc.）を用い
写像 χ, χ,̂ g を近似することで，ミニマルクエット乱流の低次元化モデル（dh ≪ d）を構成．
du∈ℝd ∼ O(105)h = χ(u) ∈ ℝdh

χ : ℝd → ℝdhχ ̂ : ℝdh → ℝdu3
u1
du= f(u)dtu2̂χ(χ(u))∈ ℝdh1h2

u3
u1u2

dh= g(h)dt19データ駆動型低次元モデル
ニューラルネット（autoencoder, Neural ODE, etc.）を用い
写像 χ, χ,̂ g を近似することで低次元化モデル（dh ≪ d）を構成．
低次元モデルで典型的な時間変動（SSP），統計量を再現

Linot and Graham, Dynamics of a data-driven low-dimensional model of turbulent
minimal Couette ﬂow, J. Fluid Mech. 973, A42 (2023).

20

データ駆動型低次元モデル
ニューラルネット（autoencoder, Neural ODE, etc.）を用い
写像 χ, χ,̂ g を近似することで低次元化モデル（dh ≪ d）を構成．
低次元モデルで典型的な時間変動（SSP），統計量を再現

低次元モデルでUPOを探索し，NSEのUPOsを新たに発見

エネルギー注入率Iと散逸率D
Linot and Graham, Dynamics of a data-driven low-dimensional model of turbulent
minimal Couette ﬂow, J. Fluid Mech. 973, A42 (2023).

21低次元化とその限界
自由度𝒪(105)のDNSが“15個の有効自由度”でモデル化可能“after dh ≳ 15 the joint PDFs are in excellent agreement with the DNS”.

エネルギー注入率Iと散逸率Dの確率密度関数:P(I, D)
Linot and Graham, Dynamics of a data-driven low-dimensional model of turbulent
minimal Couette ﬂow, J. Fluid Mech. 973, A42 (2023).

22低次元化とその限界
乱流の予測不可能性の起源である軌道不安定性，特に正のLyapunov指数 λnLE (n = 1,2,3,4)まで定量的に近似可能．
ただし，やはりdh ≳ 15は必要（低次元化の限界）．

λnLE

Lyapunov指数 λnLE (dh = 18固定)

正のLyapunov指数 λnLEのdh依存性

Linot and Graham, Dynamics of a data-driven low-dimensional model of turbulent
minimal Couette ﬂow, J. Fluid Mech. 973, A42 (2023).

23低次元化とその限界
MLによる低次元化限界はNSEのアトラクタ次元が決定（力学系的性質によって理解可能 / 限界と可能性の示唆）

λnLE

Lyapunov指数 λnLE (dh = 18固定)

正のLyapunov指数 λnLEのdh依存性

Linot and Graham, Dynamics of a data-driven low-dimensional model of turbulent
minimal Couette ﬂow, J. Fluid Mech. 973, A42 (2023).

24おしながきカオスと乱流カオスとLyapunov指数，Lyapunov次元
Blowout分岐と横断Lyapunov指数
乱流の力学系的解釈

データ駆動型乱流研究
MLによる乱流モデルの低次元化限界（Lyapunov次元）[1]
乱流の小スケール渦の隷属性とデータ同化限界（横断Lyapunov指数）[2]まとめ[1] Linot and Graham, Dynamics of a data-driven low-dimensional model of turbulent minimal Couette ﬂow, J. Fluid Mech. 973, A42 (2023).
[2] Inubushi, Saiki, Kobayashi, and Goto, Characterizing small-scale dynamics of Navier-Stokes turbulence with transverse Lyapunov exponents: A dataassimilation approach, Phys. Rev. Lett. 131, 254001 (2023).

25

データ同化(DA)
力学モデル＋観測データ⇨状態推定の高精度化
シミュレーションを用いた予測においては必須

乱流におけるDAの（軌道）安定性による理解？
不確定性の指数的増大 vs. DA

unu(0)[1][1] 辻野智紀，データ同化による気象衛星ひまわり８号風速推定データの台風強度予報へのインパクト

ϕ t(u(0))

u2u1

26

乱流の渦構造とDA

Fujino, Motoori, and Goto, JFM 975 (2023)
後藤 晋「発達した乱流―エネルギーカスケードをめぐって」（日本物理学会誌，73巻，2018）

DAの目標：乱流の小スケール構造の推定
エネルギーカスケードの素課程（cf. Goto, Saito, and Kawahara, (2017))

大スケール構造のデータを用いて，小スケール構造の推定は可能？条件/Raynolds数依存性？
27

DAの問題設定
準備E1(k)u: trueu (true) ，ũ (approx.)
分解: u = p + q, ũ = p̃ + q̃,(p :=∑|k|<kauk̂ e ik⋅x, q :=∑|k|≥kapuk̂ e ik⋅x).

uはNSEの解：
dpdq= F(p, q),= G(p, q)dt
dtqkkaE2(k)ũ: approx.

仮定p(t) (∀t ≥ 0)は使用可能（観測可能）
観測誤差・モデル化誤差なしq̃p̃kak28

DAの問題設定
目標：q(t)の推定E1(k)u: true

連続DA（直接置換法）
p̃(t) ≡ p(t) (∀t ≥ 0)

p

dp
dq
= F(p, q),= G(p, q)
( dt)
dt
先行研究

q{dq̃
= G(p, q̃)
dtkka

Obs. dataE2(k)ũ: approx.

Kolmogorov長ka > 0.2/η（十分な大スケールの情報があれば）∥u(t) − ũ(t)∥ → 0 (t → ∞)∴ ∥q(t) − q̃(t)∥ → 0 (t → ∞)DA成功（小スケールは隷属/同期）q̃p̃kak29

先行研究
Yoshida, Yamaguchi, and Kaneda, PRL (2005)
直接置換DA・非定常外力（negative viscosity）

Lalescu, Meneveau, and Eyink, PRL (2013)16η

直接置換DA・定常外力(Kolmogorov)・直方体領域

隷属/同期

Zhang, Dong, and Abdullah, JFM (2020)
４次元変分DA（4DVAR）・定常外力(Kolmogorov)

Leoni, Mazzino, and Biferale, PRX (2020)
Fourier nudging DA・ランダム非定常外力

DAのアルゴリズムや外力の詳細によらずに，kaη > 0.2であれば小スケールは隷属化．
⇨NSEの性質として理解可能？
30

周期箱乱流のDA実験

Qka := I − Pka, Pkav :=∑vk̂ e ikx

f(x, y, z) = [+cos x sin y, − sin x cos y,0]|k|<ka∂tu + (u ⋅ ∇)u = ∇π + νΔu + f, ∇ ⋅ u = 0{∂tq̃ + Qka((p + q̃) ⋅ ∇(p + q̃)) = ∇Qkaπ̃ + νΔq̃ + Qka f, ∇ ⋅ q̃ = 0
10000

E(k)/(ϵ 2/3η 5/3)

1000

9(a)

8k −5/3

76

100

5600(b)10700

100.1kη

1

1000

-10

10 -150.1

900kaη = 0.17

110 -51800

ΔE(t) := ∥u(t) − u(t)∥ = ∥q(t) − q̃(t)∥600700

800

900100031周期箱乱流のDA実験

Qka := I − Pka, Pkav :=∑vk̂ e ikx

f(x, y, z) = [+cos x sin y, − sin x cos y,0]|k|<ka∂tu + (u ⋅ ∇)u = ∇π + νΔu + f, ∇ ⋅ u = 0{∂tq̃ + Qka((p + q̃) ⋅ ∇(p + q̃)) = ∇Qkaπ̃ + νΔq̃ + Qka f, ∇ ⋅ q̃ = 0
10000

E(k)/(ϵ 2/3η 5/3)

1000

9(a)

8k −5/3

76100

5600(b)

10700

10-10

10 -150.10.1kη

1900

1000kaη = 0.17kaη = 0.20

1
10 -5

1800ΔE(t) := ∥u(t) − u(t)∥ = ∥q(t) − q̃(t)∥
600

700

8009001000

32周期箱乱流のDA実験

Qka := I − Pka, Pkav :=∑vk̂ e ikx

f(x, y, z) = [+cos x sin y, − sin x cos y,0]|k|<ka∂tu + (u ⋅ ∇)u = ∇π + νΔu + f, ∇ ⋅ u = 0{∂tq̃ + Qka((p + q̃) ⋅ ∇(p + q̃)) = ∇Qkaπ̃ + νΔq̃ + Qka f, ∇ ⋅ q̃ = 0
10000

E(k)/(ϵ 2/3η 5/3)

1000

9(a)

8k −5/3

76100

5600(b)10700

110 -51

10-10

10 -150.10.1kaη = 0.23800

1000kaη = 0.17kaη = 0.20

ΔE(t) := ∥u(t) − u(t)∥ = ∥q(t) − q̃(t)∥kaη > 0.2 ⇒ ΔE(t)600→ 0 ( 700∴ q2(t) →800q1(t))
kη
1

900

900

1000

大スケールの運動が小スケールの運動を決定（小スケール運動の隷属性）
33

NSE system∇ × u (u = p + q)

DA system∇ × ũ (ũ = p + q̃)横断Lyapunov指数を用いた特徴付け

35

同期多様体とその安定性
NSE+DA系：一方向結合力学系（skew-product）
dpdq= F(p, q),= G(p, q)dtdtdq̃= G(p, q̃)
dt(NSE, base)(DA, ﬁber)pqq̃

同期（不変）多様体：ℳ = {(p, q, q̃) | q = q̃}Fujisaka and Yamada, Prog. Theor. Phys. 69 (1983).
Ott and Sommerer, Phys. Lett. A 188 (1994)

同期多様体ℳの漸近安定性 ⟺小スケールの隷属性 ∥q(t) − q̃(t)∥ → 0 (t → ∞)36

同期多様体とその安定性
安定性は横断Lyapunov指数で判定可能∂Gidδq= [DqG(p, q)]δq [DqG]ij =変分方程式
dt∂qj )(（ℳからの無限小摂動δqをq̃ = q + δqとしてd q̃
= G(p, q̃)に代入し線形化）
dt

横断Lyapunov指数

1 ∥δq(T )∥λ⊥(ka) := lim ln
T→∞ T∥δq(0)∥(∥δq(T )∥ ∝ ∥δq(0)∥e λ⊥T )kaに依存
長時間平均量であり，λ⊥(ka) < 0なら乱流のどの初期状態からでもDAは成功．
DA（やML）のアルゴリズムによらないNSEの特徴量（変分方程式：∂tδq + Qka{u ⋅ ∇δq + δq ⋅ ∇u} = − ∇δπ + νΔδq）
37

同期多様体とその安定性
0.150.1λ⊥(ka)τ

0.05
9(a)

87065600(b)

-0.05

700

800

10 -510-10600-0.1501000ka η = 0.23

10 -15-0.1900ka η = 0.17ka η = 0.20

10.05

700

800

0.1

9000.15kaη

1000

0.2

0.250.3

乱流の小スケール（kaη > 0.2）は隷属的（NSEの性質）
横断Lyapunov指数の2つの極限，Reynolds数依存性？
39

横断Lyapunov指数 λ⊥(ka) (ka → + ∞)0.15

0.1λ⊥(ka)τ

0.05

0λ⊥(ka)τ = − (kaη)2

-0.05

-0.1

-0.15

00.05

0.1

0.15kaη

0.2

0.25

0.3

変分方程式：∂tδq + Qka{u ⋅ ∇δq + δq ⋅ ∇u} = − ∇δπ + νΔδqka → + ∞において∂tδq ≃ νΔδqとすると，λ⊥τ ∝ − (kaη)241

横断Lyapunov指数 λ⊥(ka = 0)0.15

λτ0.1λ⊥(ka)τ

0.05
0λ⊥(ka)τ = − (kaη)2

-0.05
-0.1
-0.15

00.05

0.1

0.15kaη0.2

0.25

0.3ka = 0のとき，Qka = I．λ⊥(ka = 0) = λ （最大Lyapunov指数）∵ 変分方程式：∂tδq + Qka{u ⋅ ∇δq + δq ⋅ ∇u} = − ∇δπ + νΔδq43

（再掲）乱流の力学系的解釈
乱流 ⊂ NSEによって定まる力学系のカオス
乱流のLyapunov指数： λ ∝ τ −1 (τ = (ν/ϵ)1/2: Kolmogorov time)[1]
近年の数値的研究（DNS）[2-4]：λτはわずかにRe数に依存して増大“we ﬁnd that it increases with the Reynolds number, but surprisingly faster than what ispredicted on dimensional grounds”[2][2][1] Ruelle, Phys. Lett. A 72 (1979)[2] Boﬀetta and Musacchio, PRL 119 (2017)[3] Mohan, Fitzsimmons, and Moser, PRF 2 (2017)[4] Arjun and Ho, PRL 120 (2018)44

横断Lyapunov指数のReynolds数依存性0.15

λτ0.1λ⊥(ka)τ

0.05

Re = 570

0-0.05
-0.1
-0.15

00.05

0.1

0.15kaη

0.2

0.25

0.345

横断Lyapunov指数のReynolds数依存性
λτ0.150.1λ⊥(ka)τ

0.05

Re = 570

Re = 1400

0-0.05

-0.1

-0.15

00.05

0.1

0.15

0.2

0.25

0.3kaη

最大Lyapunov指数と横断Lyapunov指数の関係：λ⊥(ka = 0) = λ
臨界スケールk*a η = 0.2は定数でなく，Reynolds数と共に増大

46まとめ
データ駆動型乱流研究の基礎付けに力学系理論は有効MLによる乱流モデルの低次元化限界（Lyapunov次元）[1]
乱流の小スケール渦の隷属性とデータ同化限界（横断Lyapunov指数）[2]

今後の展開と課題

MLによる乱流モデルの実現可能性（kη ≃ 0.2の“壁”）[3]
真にシミュレーションが必要な現象は一般に学習データが入手困難．
学習データがない領域には力学的知見を用いた“外挿”が必須．
（乱流のエネルギーカスケードに基づく転移学習[4]）データ駆動型GFD研究?（w/ Prof. Caulﬁeld@DAMTP，2024年3月〜)[1] Linot and Graham, Dynamics of a data-driven low-dimensional model of turbulent minimal Couette ﬂow, J. Fluid Mech. 973, A42 (2023).
[2] Inubushi, Saiki, Kobayashi, and Goto, Characterizing small-scale dynamics of Navier-Stokes turbulence with transverse Lyapunov
exponents: A data assimilation approach, Phys. Rev. Lett. 131, 254001 (2023).
[3] Matsumoto, Inubushi, and Goto, Feasibility of turbulence model by machine learning (in preparation).
[4] Inubushi and Goto, Transfer learning for nonlinear dynamics and its application to ﬂuid turbulence, Phys. Rev. E 102, 043301 (2020).

! inubushi@rs.tus.ac.jp " https://www.rs.tus.ac.jp/~inubushi

47