数理解析研究所講究録
第 1742 巻 2011 年 33-41

33Kolmogorov 流の乱流化と共変 Lyapunov 解析
犬伏 正信 *1, 小林 幹，
竹広 真一 ，
山田 道夫

京都大学数理解析研究所
M. INUBUSHI, M. U. KOBAYASHI, S. TAKEHIRO, M. YAMADA
Research Institute for Mathematical Sciences, Kyoto University

概要

双曲性は力学系理論における重要な概念の一つであるが，Navier-Stokes 方程式に従う流体運動のような物理系が双曲的であるか否かは明らかでない．ここでは Ginelli らによって提案された共変 Lyapunov 解析 [6] を用いて，2 次元トーラス上の Navier-Stokes 方程式(Kolmogorov 流) の双曲性を数値的に調べた．共変 Lyapunov 解析によって，解軌道に沿った安定多様体と不安定多様体の接平面 (Oseledec 分解) を求めることが可能になったので，それらのなす角度を測ることで双曲性を評価した．その結果，カオス解が現れる Reynolds 数において Kolmogorov 流は双曲的であることが分かった．さらに，Reynolds 数の増加に伴って安定多様体と不安定多様体のなす角度は次第に小さくなり，ある Reynolds 数において非双曲的となることが示唆された．この双曲性の変化は物理的な性質の変化と何らかの関係がある可能性がある．ここでは物理的な性質として，渦度の時間相関とエネルギー散逸率の時間変動の統計的性質について述べる．

1

はじめに一般に気体や液体などの流れは粘性を小さくすると (Reynolds 数を大きくすると) 外力や境界条件の与え方に依らない，普遍的な性質をもつ乱流状態になる．この流体乱流の普遍的な性質は観測や室内実験・数値実験によって支持されているが，それらに対する物理的理解が充分であるとは言えない．

乱流を理解する一つの試みとして力学系的なアプローチがある．これは乱流の物理的な性質を力学系的な視点から捉え直し，力学系の概念や解析手法を用いて乱流の理解を深めようとするものである．このアプローチに沿って，発達した乱流のシェルモデルのLyapunov スペクトル [17] や，シアー流のモデルのカオスの縁 (edge of chaos) [4], 低Reynolds 数 Navier-Stokes 乱流の不安定周期軌道 [7, 9] などが調べられ，いくつかの重要な示唆が与えられてきた．しかしながら乱流の力学系的な性質や，物理的性質との関係には不明な点が数多く残されているのが現状である．
力学系の基本的な性質の 1 つに双曲/非双曲性がある．双曲性とは，力学系の作用によって相空間が不安定な方向と安定な方向に分解されること (Oseledec 分解), すなわち安定多様体と不安定多様体が横断的に (角度をもって) 交わることをいう．力学系が双曲性を持てば理論的な取り扱いが比較的容易になる．また双曲性は構造安定性など力学系理$*1$E-mail: minubush@kurims.kyoto-u.ac.jp

34
論の中心的な概念とも密接に関係している重要な性質である．

双曲性は，解軌道上での安定/不安定多様体のなす角度を測ることで定量的に評価することが出来る．一般に安定/不安定多様体は数値的に求めることが非常に困難なため実際の計算例はほとんど無かったが，2007 年に Ginelli らが共変 Lyapunov 解析とよばれる数値アルゴリズムを開発したことにより，安定/ 不安定多様体に接するベクトル (共変Lyapunov ベクトル) を数値的に求めることが可能になった [6]. このベクトルを基に，安定/不安定多様体のなす角度の分布を計算することで，双曲性を定量的に調べることが出来る．

uptsov らは，このアルゴリズムを用いて結合 Ginzburg-Landau 方程式系の双曲/ 非双曲性を議論した [11]. 彼らは，パラメータを変化させたときの系の双曲性の程度の変化を調べ，3 番目の Lyapunov 指数が正になるパラメータと，系が非双曲的になるパラメータが一致しており，さらにそれに伴って時空カオスが現れると主張している．これは一般に双曲性の変化が物理現象に強く反映されることを示唆しており，興味深い結果である．
このアルゴリズムは従来の Lyapunov 解析に比べて大きな計算量を必要とするため，多次元偏微分方程式における実行は容易ではない．本研究では 2 次元トーラス上での非圧縮性 Navier-Stokes 方程式 (Kolmogorov 流) における乱流化過程に注目し，系の双曲/$K$非双曲性の程度が Reynolds 数とともにどのように変化するのかを調べる．また，双曲性の変化と物理的な性質の変化の関係について調べる．ここでは物理的性質として，渦度の時間相関とエネルギー散逸率の時間変動の統計的性質に着目した結果を述べる．

2問題設定と数値解析法
共変 Lyapunov 解析の方法に従って系の双曲性を数値的に評価するためには，解軌道とLyapunov ベクトルを求め，それらを用いて安定/不安定多様体のなす角度を測る必要がある．ここではまず Kolmogorov 流と呼ばれる 2 次元トーラス上 Navier-Stokes 方程式系について説明し，その数値計算法について述べる．続いて Lyapunov ベクトルの性質とその数値計算法 (共変 Lyapunov 解析) を概説する．

Kolmogorov 流は 2 次元トーラス $T^{2}((x, y)\in[0,2\pi]\cross[0,2\pi])$ 上の非圧縮 NavierStokes 方程式で記述される流れであり，渦度方程式は$\partial_{t}(+u\cdot\nabla\zeta=\frac{1}{R}(\Delta\zeta-n^{3}$$C$

O $ny)$ ,$S$(1)である．ここで $u=u(x, y, t)=(u, v)$ は速度，は渦度，は Reynolds 数， は外力の波数を表す．支配方程式 (1) は自明解と呼ばれる定常解$\zeta=-n\cos ny$ を持つ．以下，自明解の臨界 Reynolds 数 (線形不安定になる Reynolds数 を $R_{cr}(=n\sqrt{2})$と書く．外力波数が $n=1$ のとき，自明解は大域的に漸近安定で$\zeta$$=\zeta(x, y, t)=\partial_{x}v-\partial_{y}u$$R$$n$$)$あることが Iudovich によって証明されているので [8], 本研究では外力波数を $n=2$ と35する．この外力波数の下で Reynolds 数を大きくしていくと，Kolmogorov 流はカオス的になることが数値的に示されている [2](3 節参照).

Kolmogorov 流の解は，渦度をスペクトル法を用いて展開し，$\zeta(x, y, t)=\sum_{k=-K}^{K,L},$ $\zeta_{kl}(t)e^{i(kx+ly)}$,(2)

渦度方程式 (1) を直接計算することで求めた．alias 誤差は 2/3 法を用いて除去した．こ$(K=L=7)$ . 本研究で得られた結果の数値的収束は，より高い解像度の計算を行うことで確認した (確認には格子点数$30\cross 30$ ,の計算を用いた). 本研究では Kolmogorov 流の状態は実ベク
切断波数
こで格子点数は $24\cross 24$ , 切断波数は $7\cross 7$ とした$8\cross 8$

トル$(\zeta_{-K,-L}^{R}(t), \zeta_{-K,-L}^{I}(t), \cdots, \zeta_{K,L}^{R}(t) , \zeta_{K,L}^{I}(t))\in \mathbb{R}^{N}(N=2(2K+1)(2L+1))$は複素 Fourier 係数によって表されるものとする．ここで部を表す．渦度は実数値関数なので，Fourier 係数は$\zeta_{kl}^{R}(t),$$\zeta_{kl}^{I}(t)$$\zeta(z, y, t)$

す (ここである$c^{*}$

は$c\in \mathbb{C}$$\zeta_{kl}(t)$.

$( \partial_{t}\int_{T^{2}}\zeta dxdy=0)$.

さらに$\zeta_{00}$(3)の実部と虚$\zeta_{kl}(t)=\zeta_{-k-l}^{*}(t)$の複素共役を表す),を満たは Kolmogorov 流の保存量で以上の制約により本研究で扱う Kolmogorov 流の力学系としての次元 (自由度) は $N=(2K+1)(2L+1)-1=224$ である．時間積分には 4 次のRunge-Kutta 法を用$\Downarrow\rangle$, 時間刻みは $\triangle t=5.0\cross 10^{-3}$とした．初期条件は自明解に微小な摂動を加えた形として与えた．共変 Lyapunov 解析には，解がアトラクタに充分近付いたと判断した $t=1.00\cross 10^{4}$ から $t=17.0\cross 10^{4}$ までの時系列データを用いた．共変Lyapunov ベクトルは解軌道に沿って時間 $t=1$ おきに計算し，安定多様体と不安定多様体の角度を測った．

Lyapunov ベクトルとその数値計算法を概説するために， 次元離散力学系 $x_{m}=$$F(x_{m-1})$ を考える．Lyapunov ベクトル $v_{m}^{j}(j=1,2, \cdots N)$ は解軌道上の安定 (不安定)方向を指し示す接ベクトルで，そのノルムの指数的縮小 (拡大) 率が Lyapunov 指数$N$$\lambda_{j}$$(\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{N})$となるものである，すなわち$\lambda_{j}=\lim_{marrow\infty}\frac{1}{m}\log\frac{||v_{m}^{j}||}{||v_{0}^{j}||}$,$v_{m+1}^{j}=dF(x_{m})v_{m}^{j}$を満たす接ベクトルである．ここで $dF(x_{m})$ は写像 $F$ の点$x_{m}$,(4)における Jacobi 行列である．従来の Lyapunov 解析を用いれば Lyapunov 指数は計算することが出来るが，対応する Lyapunov ベクトルは直交化された接ベクトル (Gram-Schmidt ベクトル) が得られるのみで，Oseledec 分解を与える本来の定義通りの Lyapunov ベクトルは得られなかった [15]. Gram-Schmidt ベクトルは内積の定義に依存し，時間正方向/負方向で異な

36
182

$\lambda_{i}$1819202122232425$R/R_{C\ulcorner}$. The horizonFig.1 Lyapunov exponents
tal axis is Reynolds number $R/R_{cr}(18.0\leq R/R_{cr}\leq 25.0)$ . Inset; Lyapunov.
dimension$\lambda_{i}(i=1,2, \cdots 5, \lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{5})$$D_{L}$る Lyapunov 指数を与える接ベクトルである．他方，Lyapunov ベクトルは内積の定義に依らず，正方向/負方向で (符号の正/負を除いて) 同じ Lyapunov 指数を与える接ベクトルである．

系の双曲性を評価するために，Lyapunov ベクトルを用いて局所安定/不安定多様体のなす角度$\theta$$\theta=\cos^{-1}$$\max$(5)$(u, v)$ ,$u\in E^{\epsilon},v\in E^{u}||u||--||v||=1$

は正/負の Lyapunov 指数に対応する Lyapunov
ベクトルで張られる不安定/安定部分空間である．この角度は主角度 (principal angle)の方法を用いて求めた [3]. この角度 を解軌道に沿って充分長く計算することで，アト) が近似的に得られる．
ラクタの安定/不安定多様体のなす角度の分布 (確率密度関数を計算する必要がある．ここで$E^{u},$ $E^{s}$$\theta$$P(\theta)$

この確率密度関数が角度ゼロから離れていれば $(P(O)=0)$ , アトラクターは双曲的であることが分かる．また，確率密度関数が角度ゼロに値を持てば $(P(O)>0)$ , アトラクターは非双曲的になっていることがわかる．このように，Lyapunov ベクトルを用いてアトラ
クタの双曲性を評価することが出来る．

37

$0$

002

$0\alpha$$0\alpha$ou$0$$oo\epsilon$$0\alpha$$0\aleph$$0\infty$

$0l$

ou

$0\alpha$$\theta\alpha$

$0\{$

am

$0\Re$

ot

ゆ$0t$

$0$

$6\partial_{\grave{4}}$

$0$

$ooa$

ou

ゆ

ゆかさ

ゆ $\alpha$

なゆ$\theta$$0\mathfrak{l}$$\Theta$

Fig.2 Close-up $(0\leq\theta\leq 0.1[rad])$ of Probability density functions of the angle between the local stable and unstable manifolds
at Reynolds number$R/R_{cr}=20.0,21.0,22.0,23.0,24.0hom$ top to bottom (linear-log plot).
$P(\theta)$

3

共変 Lyapunov 解析
Kolmogorov 流に対して共変 Lyapunov 解析を行った結果を示す．Fig. 1 に Lyapunov

指数 $\lambda_{i}(i=1,2, \cdots 5)$ の Reynolds 数依存性を示す．Lyapunov 指数は Reynolds 数と$R/R_{cr}=18.0$ では準周期解が安定であるが共に単調に増大することがわかる．特に，$R/R_{cr}\simeq 18.2$ で不安定化し，カオス解が現れること$(\lambda_{1}=\lambda_{2}=0, \lambda_{i}\leq 0for^{\forall}i),$$(\lambda_{1}>0)$,さらに $R/R_{cr}>23.0\sim$ では 2 番目の正の Lyapunov 指数が現れることがわかる $(\lambda_{1}>\lambda_{2}>0)$ .

Fig. 1 の内挿図は Lyapunov 次元$K= \max\{m|\sum_{j=1}^{m}\lambda_{j}\geq 0\})$$D_{L}=K+ \frac{1}{|\lambda_{K+1}|}\sum_{i=1}^{K}\lambda_{i}$(ここでである．Lyapunov 次元は $3.5<D_{L}<\sim\sim 5.5$ であることがわかる．本研究では Reynolds 数領域 20. $0\leq R/R_{cr}\leq 24.0$ において共変 Lyapunov 解析を行った．

Kolmogorov 流の双曲性を評価するために，解軌道に沿って測った局所安定/不安定多様体のなす角度 の確率密度関数を示す (Fig.2). Reynolds 数は上から順に$\theta$$R/R_{cr}=20.0,21.0,22.0,23.0,24.0$$P(\theta)$てある．双曲性の議論には角度ゼロの確率 $P(O)$ が重38$0$

500

1000
$\tau$. EnlargeFig.3 The ensemble averaged time-correlation functions
ment of the time-correlation functions in $0\leq\tau\leq 10$ is shown in the inset.
$\langle\rho(\tau)\rangle/\langle\rho(0)\rangle$[rad]
要なので，Fig.2 には得られた確率密度関数の角度ゼロ付近の拡大図を示した．Reynolds 数が小さいとき $(R/R_{cr}\simeq 20.0)$ は，局所安定/不安定多様体のなす角度はゼロから離れていることがわかる $(P(O)=0)$ . つまりアトラクタは双曲的である．しかし，Reynolds 数が大きくなるにつれて次第に局所安定/不安定多様体のなす角度は小さな値をとるようになる．さらにある Reynolds 数 $(R/R_{cr}\simeq 23.0)$ で確率密度関$(0\leq\theta\leq 0.1$$)$

数は角度ゼロに確率を持つようになる $(P(O)>0)$ , これはアトラクタが非双曲的になっていることを示唆する．また，双曲性が破れる Reynolds 数 $(R/R_{cr}\simeq 23.0)$ は 2 番目の正の Lyapunov 指数が現れる Reynolds 数に近いことがわかる．先述の通り，Kuptsov
らは結合 Ginzburg-Landau 方程式系において双曲性が破れるパラメタと 3 番目の正のLyapunov 指数が現れるパラメタが一致していると主張している [11]. index が異なるものの，これらの結果は系の不安定化と双曲性の破れに関係があることを示唆しているという点で類似した結果といえる．

4

物理的性質との関係
前節では Kolmogorov 流の双曲性が Reynolds 数の増大に伴って変化することがわかった．本節では，渦度の時間相関とエネルギー散逸率の時間変動に着目し，双曲性の変化と物理的な性質の変化との関係を調べる．

39

Fig.4 Joint probability density functions of the energy dissipation rate and theangle . The horizontal axis is the angle , the vertical axis is the energy dissipation rate and the contour is the joint PDF at $R/R_{cr}=20.0(a),$ $22.0(b),$ $24.0(c)$ .
$\epsilon$$\theta$$\theta$$\epsilon$

まず $(x’, y’)=(\pi/4, \pi/4)$ における渦度 $\zeta(t)=\zeta(x’, y’, t)$ の時間相関$\rho(\tau)$.

(6)$\rho(\tau)=\lim_{Tarrow\infty}\frac{1}{T}\int_{0}^{T}\hat{\zeta}(t)\hat{\zeta}(t-\tau)dt$

について述べる．Fig.3 に$\langle\rho(\tau)\rangle/\langle\rho(0)\rangle$を示す．ここで$\hat{\zeta}(t)$は平均値からの変動$\hat{\zeta}(t)=$を表し， は異なる 60 個の初期値を用いたアンサンブル平均を表す．充分長い平均時間として $T=2.0\cross 10^{5}$ を用いた．Reynolds 数 $R/R_{cr}=20.0,21.0,22.0,23.0,24.0$$\zeta(t)-\overline{\zeta}$$\langle\cdot\rangle$における時間相関関数を示してある．Fi 3 の内挿図は相関関数の短い時間 $(0\leq\tau\leq 10)$の拡大図である．これより，短い時間の相関は Reynolds 数の増大に伴って徐々に失われることがわかる．他方，Fig.3 に示された長時間相関 $(\tau\sim>500)$ の値に見られるように，長時間相関は 2 つのグループ $(R/R_{cr}=20.0,21.0$ と $R/R_{cr}=22.0,23.0,24.0)$ にわかれる．長時間相関に変化が起きる Reynolds 数 $(R/R_{cr}\simeq 22.0)$ は，アトラクタの双曲性が破れる Reynolds 数 $(R/R_{cr}\simeq 23.0)$ に近く，これらの変化の間には何らかの関係があることが示唆される．

最後に，エネルギー散逸率 と局所安定/不安定多様体のなす角度 の関係について調べた．エネルギー散逸率とその時間空間的なゆらぎは，Kolmogorov の次元解析による乱$\theta$$\epsilon$

流理論，及び高次統計量に見られる理論からのずれ (間欠性) を考える上で重要な物理量である．エネルギー散逸率 と局所安定/不安定多様体のなす角度 の結合確率密度関数を Fig.4 に示す．Reynolds 数はそれぞれ (a) $R/R_{cr}=20.0,(b)R/R_{cr}=22.0,(c)$$R/R_{cr}=24.0$ である．Reynolds 数 $R/R_{cr}=24.0$ における結合確率密度関数
か$\theta$$\epsilon$$P(\theta, \epsilon)$$P(\theta, \epsilon)$

ら，エネルギー散逸が強く起こっているときは角度$\theta$は小さい，という関係があることがわかる．この関係は $R/R_{cr}=23.0,25.0$ でも同様に見られた．この結果により．系の双曲性が破れているときはエネルギー散逸率$\epsilon$と局所安定/不安定多様体のなす角度$\theta$に上に述べたような関係があることがわかった．同様の関係は，乱流のシェルモデル (GOYモデル) においても見られる [10]. GOY モデルは発達した 3 次元乱流のモデルであり，

40
力学系的性質 (アトラクター次元$D_{L}$など) を見ても本研究で扱った Kolmogorov 流とは大きく異なる．それにも関わらず，共通の性質を示すことは興味深い．この関係の物理的な意味を明かにすること，他の一般の散逸力学系で同様の関係が成り立っているかを調べることは今後の課題である．
謝辞: 本研究会で講演の機会を与えて下さいました新居俊作先生に深く感謝致します．

参考文献[1] Z. Arai, On Hyperbolic Plateaus of the Henon Map, Experimental Mathematics
16, 2, pp.181-188 (2007).
[2] D. Armbruster, B. Nicolaenko, N. Smaoui, and P. Chossat, Symmetries anddynamics for 2-D Navier-Stokes flow, Physica $D95$ , pp.81-93 (1996).
[3] A. Bjorck, and G. H. Golub, Numerical methods for computing angles between
linear subspages, Mathematics of computation, vol. 27,Number 123, pp.579-594,174101 (4) (1973).
[4] J. D. Skufca, J. A. Yorke, and B. Eckhardt, Edge of Chaos in a Paralle117 Shear
Flow, Phys. Rev. Lett. 96, 174101(4) (2006).
[5] J. P. Eckmann, and D. Ruelle, Ergodic theory of chaos and strange attractors,
Reviews of Modern Physics, vol.57,No.3, pp.617-656 (1985).
[6] F. Ginelli, P. Poggi, A. Turchi, H. Chate, R. Livi, and A. Politi, Characterizing Dynamics with Covariant Lyapunov Vectors, Phys. Rev. Lett. 99,130601(4)(2007).
[7] T. Itano, and S. Toh, The Dynamics of Bursting Process in Wall Turbulence,
Journal of Physical Society of Japan, Vo170, No.3, pp.703-716 (2000).
[8] V. I. Iudovich, Example of the generation of a secondary stationary or periodic
flow when there is loss of stability of the laminar flow of a viscous incompressible
fiuid, J. Appl. Math. Mech. 29, pp.527-544 (1965).
[9] G. Kawahara, and S. Kida, Periodic motion embedded in plane Couette turbulence: regeneration cycle and burst, J. Fluid. Mech. 449, pp.291-330 (2001).
[10] M. U. Kobayashi and M. Yamada, in preparation.
[11] P. V. Kuptsov, and S. P. Kuznetsov, Violation of hyperbolicity in a diffusive
medium with local hyperbolic attractor, phys. rev. $E80$ , 016205(11) (2009).
[12] L. D. Meshlkin, and Y. G. Sinai, Investigation of the stability of a stationary
solution of a system of equations for the plane movement of an incompressible
viscous liquid, Journal of applied Mathematics. 25, pp.1700-170518 (1962).
[13] N. Platt, L. Sirovich, and N. Fitzmaurice, An investigation of chaotic Kolmogorov
flows, Physics of Fluids A,Vol. 3, No. 4, pp.681-696 (1990).
[14] Y. Saiki and M. U. Kobayashi, Numerical Identffication of Nonhyperbolicity of
the Lorenz System through Lyapunov Vectors, JSIAM Letters Vol.2, pp.107-110(2010)..
[15] I. Shimada, and T. Nagashima, A numerical approach to ergodic problem of dissipative dynamical systems, Progress of theoretical physics,vol.61,No.6, pp.16051616 (1978).
[16] S. Takehiro, M. Odaka, K. Ishioka, M. Ishiwatari, Y. Hayashi, SPMODEL: A41
Series of Hierarchical Spectral Models for Geophyiscal Fluid Dynamics, Nagare
Multimedia, http: $//www$ .nagare. or. $jp/mm/2006/spmodel/(2006)$ .

[17] M. Yamada, and K. Ohkitani,Asymptotic formulas for the Lyapunov spectrum offully developed shell model turbulence, Phys. Rev. E. 57, p.p.6257-6260 (1998).