% 表題 回転球面上での 2 次元順圧流体の帯状流安定性問題の定式化 % % 履歴 2005/07/21 竹広 真一 新規作成 % 2005/08/11 竹広 真一 Rayleigh's criterion 追加 % 2007/11/25 竹広 真一 参照文書追加 % % \documentclass[12pt,a4j,notitlepage]{jarticle} \usepackage{Dennou6} \usepackage{psfrag} \usepackage{ascmac} \usepackage{times} \Dtitle[回転球面順圧流体〜帯状流安定性] {回転球面上での 2 次元順圧流体\\〜帯状流安定性問題の定式化} \Dauthor[竹広 真一]{竹広 真一, SPMODEL 開発グループ} \Ddate{平成 19 年 11 月 25 日} \Dparskip \pagestyle{Dheadings} \begin{document} \maketitle \section{はじめに} この文章では, 回転球面上の 2 次元順圧系での帯状流の安定性を 数値的な時間積分により調べるための定式化と問題設定について述べる. \section{支配方程式} 支配方程式は, 回転球面上での渦度方程式である (「2 次元回転球面上での非発散順圧流体の定式化」参照): % \begin{equation} \Deqlab{渦度方程式} \DP{\zeta}{t} + \Dinv{a^2}J(\psi,\zeta) + \frac{2\Omega}{a^2}\DP{\psi}{\lambda} = (-1)^{p+1}\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p \zeta, \end{equation} % ここで $\zeta$ は渦度の動径成分, $t$ は時間, $\psi$ は流線関数, $a$ は 球の半径, $\Omega$ は球の回転角速度, $\lambda$ は経度, $\nu_p$ は高階粘性の係数である\footnote{$p=1$ で通常の粘性項になる}. $\Ddsty J(f,g)\equiv \DP{f}{\lambda}\DP{g}{\mu}-\DP{f}{\mu}\DP{g}{\lambda} = \Dinv{\cos\varphi}\DP{f}{\lambda}\DP{g}{\varphi} -\DP{f}{\varphi}\Dinv{\cos\varphi}\DP{g}{\lambda}$ はヤコビアンであり, $\varphi$ は緯度, $\mu=\sin\varphi$ は sin 緯度 である. 流線関数と速度の経度・緯度成分 $u_\lambda,u_\varphi$ との関係は % \begin{equation} \Deqlab{速度流線関数関係} u_\lambda = - \Dinv{a}\DP{\psi}{\varphi}, \quad u_\varphi = \Dinv{a\cos\varphi}\DP{\psi}{\lambda}, \end{equation} % であり, したがって渦度と流線関数の関係は % \begin{eqnarray} \Deqlab{渦度流線関数関係} \zeta &=& (\Drot\Dvect{u})_r = \Dinv{a\cos\varphi}\DP{u_\varphi}{\lambda} - \Dinv{a\cos\varphi}\DP{(u_\lambda\cos\varphi)}{\varphi} \nonumber\\ &=& \Dinv{a^2\cos^2\varphi}\DP[2]{\psi}{\lambda} + \Dinv{a^2\cos\varphi}\DP{}{\varphi} \left(\cos\varphi\DP{\psi}{\varphi}\right) \nonumber\\ &=& \Dlapla \psi. \end{eqnarray} % ここで $\Ddsty \Dlapla \equiv \Dinv{a^2\cos^2\varphi}\DP[2]{}{\lambda} + \Dinv{a^2\cos\varphi}\DP{}{\varphi} \left(\cos\varphi\DP{}{\varphi}\right)$ は半径 $a$ の球面上での 2 次元水平ラプラシアンである. \section{数値計算 : 球面調和函数展開と時間積分} \subsection{球面調和函数展開} スペクトル法による数値計算を行うために, 各物理量を球面調和函数 $Y_n^m(\lambda,\varphi)$で 展開する. すなわち, % \begin{eqnarray} \zeta(\lambda,\varphi) &=&\sum_{n=0}^\infty\sum_{m=-n}^n \tilde{\zeta}_n^m(t)Y_n^m(\lambda,\varphi), \\ \psi(\lambda,\varphi) &=&\sum_{n=0}^\infty\sum_{m=-n}^n \tilde{\psi}_n^m(t)Y_n^m(\lambda,\varphi) \end{eqnarray} % ここで, 球面調和函数は数値計算のため実関数として つぎのように定義されたものを扱う. % \begin{eqnarray} Y_n^m(\lambda,\varphi) &=& P_n^{m}(\cos\varphi)\cos (m\lambda) \quad (m \ge 0),\\ &=& P_n^{|m|}(\cos\varphi)\sin (|m|\lambda) \quad (m < 0), \end{eqnarray} % これを\Deqref{渦度方程式}に代入し, $Y_{n'}^{m'}(\lambda,\varphi)$ をかけて全球面で積分することにより, 球面調和函数の直交関係から以下の式が得られる. % \begin{displaymath} \DD{\tilde{\zeta}_n^m}{t} = - \Dinv{a^2}[J(\psi,\zeta)]_n^m - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m + \left[(-1)^{p+1}\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p \zeta \right]_n^m \end{displaymath} % ただし $[\ldots]_n^m$ は球面調和函数の $n,m$ 成分であることを表している. % 散逸項はラプラシアンが球面調和函数の固有演算子であることから \begin{eqnarray*} \left[(-1)^{p+1}\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p \zeta(\lambda,\varphi,t) \right]_n^m &=& (-1)^{p+1}\nu_{2p}\left[-n(n+1) + \frac{2}{a^2}\right]^p \tilde{\zeta}_n^m(t)\\ &=& -\nu^*(n,p)\tilde{\zeta}_n^m(t), \end{eqnarray*} % となる. ここで, 波数毎の超粘性係数 $\nu^*(n,p)$ を \begin{equation} \nu^*(n,p)\equiv (-1)^p\nu_{2p} \left[-n(n+1) + \frac{2}{a^2}\right]^p \end{equation} % と定義した($\nu^*(n,p)$ は 0 以上であることに注意). したがって, \begin{equation} \Deqlab{スペクトル渦度方程式} \DD{\tilde{\zeta}_n^m}{t} = - \Dinv{a^2}[J(\psi,\zeta)]_n^m - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m - \nu^*(n,p)\tilde{\zeta}_n^m \end{equation} \subsection{時間積分} \subsubsection{euler スキーム(1 次精度)} euler スキームの場合は時間変化項をすべて現在時間で評価して, % \begin{eqnarray*} \frac{\tilde{\zeta}(t+\Delta t)_n^m-\tilde{\zeta}(t)_n^m}{\Delta t} &=& - \Dinv{a^2}[J(\psi(\lambda,\varphi,t),\zeta(\lambda,\varphi,t))]_n^m - \frac{2\Omega}{a^2}\left[\DP{\psi(\lambda,\varphi,t)}{\lambda}\right]_n^m \\ && - \nu^*(n,p)\tilde{\zeta}_n^m(t) \end{eqnarray*} % すなわち, % \begin{eqnarray} \tilde{\zeta}(t+\Delta t)_n^m =\tilde{\zeta}(t)_n^m &+&\Delta t \times \left\{ - \Dinv{a^2}[J(\psi(\lambda,\varphi,t),\zeta(\lambda,\varphi,t))]_n^m - \frac{2\Omega}{a^2}\left[\DP{\psi(\lambda,\varphi,t)}{\lambda}\right]_n^m \right. \nonumber\\ & &\left. \hspace{20mm} - \nu^*(n,p)\tilde{\zeta}_n^m(t) \right\}. \end{eqnarray} \subsubsection{Crank-Nicolson and Adams-Bashforth スキーム(2 次精度)} 散逸項に Clank-Nicolson スキーム, それ以外の項に 2 次の Adams-Bashforth スキームを適用する場合には, 時間積分を 2 段階に分けて % \begin{eqnarray*} \frac{\tilde{\zeta}^*(t+\Delta t)_n^m-\tilde{\zeta}(t)_n^m}{\Delta t} &=& \frac{3}{2}\left[\left(\DD{\zeta}{t}\right)_{nd}(t)\right]_n^m -\frac{1}{2}\left[\left(\DD{\zeta}{t}\right)_{nd}(t-\Delta t)\right]_n^m,\\ \frac{\tilde{\zeta}(t+\Delta t)_n^m-\tilde{\zeta}^*(t+\Delta t)_n^m}{\Delta t} & & = \frac{1}{2}\left\{ \left[\left(\DD{\zeta}{t}\right)_{d}(t+\Delta t)\right]_n^m +\left[\left(\DD{\zeta}{t}\right)_{d}(t)\right]_n^m \right\}, \\ \end{eqnarray*} % ただし, $\Ddsty\left(\DD{\zeta}{t}\right)_{nd}, \left(\DD{\zeta}{t}\right)_{d}$ はそれぞれ時間変化の 非散逸項と散逸項を表している: % \begin{eqnarray*} \left[\left(\DD{\zeta}{t}\right)_{nd}(t)\right]_n^m &\equiv& - \Dinv{a^2}[J(\psi(\lambda,\varphi,t),\zeta(\lambda,\varphi,t))]_n^m - \frac{2\Omega}{a^2}\left[\DP{\psi(\lambda,\varphi,t)}{\lambda}\right]_n^m,\\ \left[\left(\DD{\zeta}{t}\right)_{d}(t)\right]_n^m &\equiv& - \nu^*(n,p)\tilde{\zeta}_n^m(t) \end{eqnarray*} % すると, % \begin{eqnarray*} \frac{\tilde{\zeta}^*(t+\Delta t)_n^m-\tilde{\zeta}(t)_n^m}{\Delta t} &=& \frac{3}{2}\left[\left(\DD{\zeta}{t}\right)_{nd}(t)\right]_n^m -\frac{1}{2}\left[\left(\DD{\zeta}{t}\right)_{nd}(t-\Delta t)\right]_n^m\\ \frac{\tilde{\zeta}(t+\Delta t)_n^m-\tilde{\zeta}^*(t+\Delta t)_n^m}{\Delta t} &=& \frac{-\nu^*(n,p)\tilde{\zeta}_n^m(t+\Delta t) -\nu^*(n,p)\tilde{\zeta}_n^m(t)}{2}, \end{eqnarray*} % したがって % \begin{eqnarray} \tilde{\zeta}^*(t+\Delta t)_n^m &=& \tilde{\zeta}(t)_n^m + \Delta t \times \left\{ \frac{3}{2}\left[\left(\DD{\zeta}{t}\right)_{nd}(t)\right]_n^m -\frac{1}{2}\left[\left(\DD{\zeta}{t}\right)_{nd}(t-\Delta t) \right]_n^m \right\}, \\ % % \tilde{\zeta}(t+\Delta t)_n^m &=& \frac{1 - \nu^*(n,p)\Delta t/2}{1 + \nu^*(n,p)\Delta t/ 2} \times\tilde{\zeta}^*(t+\Delta t)_n^m. \end{eqnarray} \subsubsection{演算子分割処理法(粘性項)と 4 次精度 Runge-Kutta スキーム} 変数変換することによって粘性項を解析的に評価し, その他の項を高次の数値積分スキームを適用する 演算子分割法をとる場合を以下に示す. 渦度のスペクトル成分について % \begin{equation} \tilde{\zeta}(t)_n^m = \hat{\zeta}(t)_n^m e^{-\nu^*(n,p)t} \end{equation} % と変数変換すれば, % \begin{displaymath} \DD{\tilde{\zeta}_n^m}{t} = \DD{\hat{\zeta}(t)_n^m}{t} e^{-\nu^*(n,p)t} - \nu^*(n,p)\hat{\zeta}(t)_n^m e^{\nu^*(n,p)t} = \DD{\hat{\zeta}(t)_n^m}{t} e^{-\nu^*(n,p)t} - \nu^*(n,p)\tilde{\zeta}(t)_n^m \end{displaymath} % であるから, \Deqref{スペクトル渦度方程式}は \begin{equation} \DD{\hat{\zeta}(t)_n^m}{t} e^{-\nu^*(n,p)t} = - \Dinv{a^2}[J(\psi,\zeta)]_n^m - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m. \end{equation} % \begin{equation} \DD{\hat{\zeta}(t)_n^m}{t} = e^{\nu^*(n,p)t}\left[ - \Dinv{a^2}[J(\psi,\zeta)]_n^m - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m \right]. \end{equation} % この左辺の時間積分を Runge-Kutta スキームで計算する. すなわち, \begin{eqnarray*} \tilde{\zeta}(t+\Delta t)_n^m &=& \hat{\zeta}(t+\Delta t)_n^m e^{-\nu^*(n,p)\Delta t},\\ \hat{\zeta}(t+\Delta t)_n^m &=& \tilde{\zeta}(t)_n^m + \Delta t \times \left( \frac{k_1}{6} +\frac{k_2}{3} +\frac{k_3}{3} +\frac{k_4}{6} \right), \\ k_1 &=& f[\zeta(t),\psi(t),0]_n^m,\\ k_2 &=& f[\zeta(t_1),\psi(t_1),\Delta t/2]_n^m\\ k_3 &=& f[\zeta(t_2),\psi(t_2),\Delta t/2]_n^m\\ k_4 &=& f[\zeta(t_3),\psi(t_3),\Delta t]_n^m\\ f[\zeta,\psi,t]_n^m &=& e^{\nu^*(n,p)t}\left[ - \Dinv{a^2}[J(\psi,\zeta)]_n^m - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m \right],\\ \tilde{\zeta}_n^m(t_1) &=& e^{-\nu^*(n,p)\Delta t/2} [\tilde{\zeta}(t)_n^m + k_1 \Delta t/2],\\ \tilde{\zeta}_n^m(t_2) &=& e^{-\nu^*(n,p)\Delta t/2} [\tilde{\zeta}(t)_n^m + k_2 \Delta t/2],\\ \tilde{\zeta}_n^m(t_3) &=& e^{-\nu^*(n,p)\Delta t} [\tilde{\zeta}(t)_n^m + k_3 \Delta t]. \end{eqnarray*} % この定式化において $\tilde{\zeta}(t)_n^m = \hat{\zeta}(t)_n^m$ であることを 用いていることに注意されたい. \subsubsection{演算子分割処理法(線形項)と 4 次精度 Runge-Kutta スキーム} 粘性項だけでなく, 全ての線形項を演算子分割処理方によって解析的に扱い, 残りの非線形項のみ数値積分して解くための定式化を記す. 今, $m$ を正の整数として, \Deqref{スペクトル渦度方程式}の $(n,m)$ 成分と $(n,-m)$ 成分の式を書き下すと % \begin{eqnarray*} \DD{\tilde{\zeta}_n^m}{t} &=& - \Dinv{a^2}[J(\psi,\zeta)]_n^m - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m - \nu^*(n,p)\tilde{\zeta}_n^m, \\ \DD{\tilde{\zeta}_n^{-m}}{t} &=& - \Dinv{a^2}[J(\psi,\zeta)]_n^{-m} - \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^{-m} - \nu^*(n,p)\tilde{\zeta}_n^{-m} \end{eqnarray*} % ここで $\Dlapla\psi=\zeta$ より $\tilde{\psi}_n^m = -\tilde{\zeta}_n^m/n(n+1)$ であることを 用いると, \begin{displaymath} \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^m = \frac{2\Omega m}{a^2 n(n+1)} \tilde{\zeta}_n^{-m},\quad \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^{-m} = -\frac{2\Omega m}{a^2 n(n+1)} \tilde{\zeta}_n^m, \end{displaymath} % となるので % \begin{eqnarray*} \DD{\tilde{\zeta}_n^m}{t} &=& - \Dinv{a^2}[J(\psi,\zeta)]_n^m - \frac{2\Omega m}{a^2 n(n+1)} \tilde{\zeta}_n^{-m} - \nu^*(n,p)\tilde{\zeta}_n^m, \\ \DD{\tilde{\zeta}_n^{-m}}{t} &=& - \Dinv{a^2}[J(\psi,\zeta)]_n^{-m} + \frac{2\Omega m}{a^2 n(n+1)} \tilde{\zeta}_n^m - \nu^*(n,p)\tilde{\zeta}_n^{-m} \end{eqnarray*} % ロスビー波の位相速度 $\Ddsty\omega_R(n,m)\equiv - \frac{2\Omega |m|}{a^2 n(n+1)}$ を用いて書き直すと, % \begin{eqnarray*} \DD{\tilde{\zeta}_n^m}{t} &=& - \Dinv{a^2}[J(\psi,\zeta)]_n^m + \omega_R(n,m) \tilde{\zeta}_n^{-m} - \nu^*(n,p)\tilde{\zeta}_n^m, \\ \DD{\tilde{\zeta}_n^{-m}}{t} &=& - \Dinv{a^2}[J(\psi,\zeta)]_n^{-m} - \omega_R(n,m) \tilde{\zeta}_n^m - \nu^*(n,p)\tilde{\zeta}_n^{-m} \end{eqnarray*} % 解きやすくするために $\tilde{\zeta}_n^m$ と $\tilde{\zeta}_n^{-m}$ を複素変数 $Z_n^m$ でまとめて扱う. すなわち $ Z_n^m \equiv \tilde{\zeta}_n^m + i \tilde{\zeta}_n^{-m} $ と置き換えると, \begin{displaymath} \DD{Z_n^m}{t} = G(\psi,\zeta)_n^m -i \omega_R(n,m) Z_n^m - \nu^*(n,p)Z_n^m, \end{displaymath} % ただし $G(\psi,\zeta)_n^m$ は非線形項を複素表示したものであり, \\ $\Ddsty G(\psi,\zeta)_n^m \equiv -\left[ \Dinv{a^2}[J(\psi,\zeta)]_n^m +i\Dinv{a^2}[J(\psi,\zeta)]_n^{-m}\right]$ である. 線形項を消去すべく $Z_n^m(t) = \hat{Z}_n^m(t) e^{-[i\omega_R(n,m)+\nu^*(n,p)]t}$ と変換すれば, \begin{eqnarray*} \DD{\hat{Z}_n^m}{t} e^{-[i\omega_R(n,m)+\nu^*(n,p)]t} &-&[i\omega_R(n,m)+\nu^*(n,p)]\hat{Z}_n^m(t) e^{-[i\omega_R(n,m)+\nu^*(n,p)]t}\\ &=& G(\psi,\zeta)_n^m -[i\omega_R(n,m) Z_n^m + \nu^*(n,p)]Z_n^m, \end{eqnarray*} したがって \begin{displaymath} \DD{\hat{Z}_n^m}{t} = e^{[i\omega_R(n,m)+\nu^*(n,p)]t} G(\psi,\zeta)_n^m \end{displaymath} % 実部と虚部に分けて $Re[\hat{Z}_n^m]=\xi_n^m, Im[\hat{Z}_n^m]=\xi_n^{-m}$ とおくと % \begin{eqnarray} \DD{\xi_n^m}{t} &=& e^{\nu^*(n,p)]t}\cos(\omega_R(n,m) t) Re[G(\psi,\zeta)_n^m] - e^{\nu^*(n,p)]t}\sin(\omega_R(n,m) t)Im[G(\psi,\zeta)_n^m]\nonumber\\ &=& e^{\nu^*(n,p)]t}\cos(\omega_R(n,m) t) \left[-\Dinv{a^2}J(\psi,\zeta)\right]_n^m - e^{\nu^*(n,p)]t}\sin(\omega_R(n,m) t) \left[-\Dinv{a^2}J(\psi,\zeta)\right]_n^{-m},\\ \DD{\xi_n^{-m}}{t} &=& e^{\nu^*(n,p)]t}\cos(\omega_R(n,m) t) Im[G(\psi,\zeta)_n^m] + e^{\nu^*(n,p)]t}\sin(\omega_R(n,m) t) Re[G(\psi,\zeta)_n^m]\nonumber\\ &=& e^{\nu^*(n,p)]t}\cos(\omega_R(n,m) t) \left[-\Dinv{a^2}J(\psi,\zeta)\right]_n^{-m} + e^{\nu^*(n,p)]t}\sin(\omega_R(n,m) t) \left[-\Dinv{a^2}J(\psi,\zeta)\right]_n^m,\\ \end{eqnarray} % 時間積分のための $\xi_n^m$, $\xi_n^{-m}$ の各積分ステップでの初期値は $Z_n^m$ と $\hat{Z}_n^m$ の関形において $t=0$ を適用すればよく, $\xi_n^m = \tilde{\zeta}_n^m$, $\xi_n^{-m} = \tilde{\zeta}_n^{-m}$ となる. 上の式を用いて数値積分し, あたらしい時間の $\xi_n^m$, $\xi_n^{-m}$ が求まれば, 対応して $\tilde{\zeta}_n^m$ がつぎのように求められる. \begin{eqnarray} \tilde{\zeta}_n^m &=& Re[Z_n^m] = Re[\hat{Z}_n^m e^{-[i\omega_R(n,m)+\nu^*(n,p)]t}]\nonumber\\ &=& Re[\hat{Z}_n^m] e^{-\nu^*(n,p)t}\cos(\omega_R(n,m) t) + Im[\hat{Z}_n^m] e^{-\nu^*(n,p)t}\sin(\omega_R(n,m) t)\nonumber\\ &=& \xi_n^m e^{-\nu^*(n,p)t}\cos(\omega_R(n,m) t) + \xi_n^{-m} e^{-\nu^*(n,p)t}\sin(\omega_R(n,m) t),\\ \tilde{\zeta}_n^{-m} &=& Im[Z_n^m] = Im[\hat{Z}_n^m e^{-[i\omega_R(n,m)+\nu^*(n,p)]t}]\nonumber\\ &=& Im[\hat{Z}_n^m] e^{-\nu^*(n,p)t}\cos(\omega_R(n,m) t) - Re[\hat{Z}_n^m] e^{-\nu^*(n,p)t}\sin(\omega_R(n,m) t)\nonumber\\ &=& \xi_n^{-m} e^{-\nu^*(n,p)t}\cos(\omega_R(n,m) t) - \xi_n^m e^{-\nu^*(n,p)t}\sin(\omega_R(n,m) t). \end{eqnarray} \section{帯状流の安定性の理論的考察} 以下では, 帯状流の安定性を理論的に考察するために 基本場(帯状流)と擾乱に分けた方程式を展開し, 帯状流が不安定になるための必要条件を導く. \subsection{非粘性の支配方程式と保存量} 理論的考察においては, 数値計算の安定のために \Deqref{渦度方程式}の右辺に導入していた高階粘性を無視した式, % \begin{equation} \Deqlab{非粘性渦度方程式} \DP{\zeta}{t} + \Dinv{a^2}J(\psi,\zeta) + \frac{2\Omega}{a^2}\DP{\psi}{\lambda} = 0, \end{equation} % を扱うことにする. 非粘性の下ではいくつかの保存量が存在する. 中でも重要なものはポテンシャル渦度保存則である. \Deqref{非粘性渦度方程式}の左辺第 3 項を変形して \begin{eqnarray*} \frac{2\Omega}{a^2}\DP{\psi}{\lambda} &=& \Dinv{a^2\cos\varphi}\DP{(2\Omega\sin\varphi)}{\varphi} \DP{\psi}{\lambda}\\ &=& \Dinv{a^2\cos\varphi}\DP{(2\Omega\sin\varphi)}{\varphi} \DP{\psi}{\lambda} - \Dinv{a^2\cos\varphi}\DP{(2\Omega\sin\varphi)}{\lambda} \DP{\psi}{\varphi}\\ &=& -\Dinv{a^2}J(2\Omega\sin\varphi,\psi) = \Dinv{a^2}J(\psi,2\Omega\sin\varphi), \end{eqnarray*} % したがって % \begin{displaymath} \DP{\zeta}{t} + \Dinv{a^2}J(\psi,\zeta+2\Omega\sin\varphi)= 0,\\ \end{displaymath} % $2\Omega\sin\varphi$ は時間によらないので, % \begin{displaymath} \DP{(\zeta+2\Omega\sin\varphi)}{t} + \Dinv{a^2}J(\psi,\zeta+2\Omega\sin\varphi)= 0, \end{displaymath} % ここにポテンシャル渦度 $q(\lambda,\varphi,t)$ を % \begin{equation} q=2\Omega\sin\varphi + \zeta \end{equation} % と定義すれば, その変化を表す式は % \begin{equation} \Deqlab{ポテンシャル渦度保存則} \DP{q}{t} + \Dinv{a^2}J(\psi,q)= 0, \end{equation} % あるいは物質微分を用いて % \begin{equation} \frac{D q}{D t } =0, \end{equation} % と表される. すなわちポテンシャル渦度は流体粒子とともにその値を 変えることなく流れていく保存量である. 2 次元順圧系ではポテンシャル渦度は慣性系からみたときの 渦度である絶対渦度に等しい. これに応じて $2\Omega\sin\varphi$ を惑星渦度, $\zeta$ を相対渦度とよぶ. 他の保存量として, エンストロフィー(渦度の2 乗)が存在する. ポテンシャル渦度保存則に $q$ をかけると % \begin{equation} \Deqlab{絶対エンストロフィー保存則} \DP{}{t}\left(\Dinv{2}q^2\right) + \Dinv{a^2}J\left(\psi,\Dinv{2}q^2\right)= 0, \end{equation} % 相対渦度に対するエンストロフィー $\zeta^2/2$ で書くならば \Deqref{非粘性渦度方程式}に $\zeta$ をかけて % \begin{displaymath} \DP{}{t}\left(\Dinv{2}\zeta^2\right) + \Dinv{a^2}J\left(\psi,\Dinv{2}\zeta^2\right) + \frac{2\Omega}{a^2}\zeta\DP{\psi}{\lambda} = 0, \end{displaymath} % 第 2 項目は \begin{eqnarray*} \Dinv{a^2}J\left(\psi,\Dinv{2}\zeta^2\right) &=& \Dinv{a^2\cos\varphi}\left[ \DP{\psi}{\lambda}\DP{}{\varphi}\left(\Dinv{2}\zeta^2\right) -\DP{\psi}{\varphi}\DP{}{\lambda}\left(\Dinv{2}\zeta^2\right) \right]\\ &=& \Dinv{a^2\cos\varphi}\DP{}{\varphi} \left[\DP{\psi}{\lambda} \left(\Dinv{2}\zeta^2\right)\right] -\Dinv{a^2\cos\varphi} \DP{}{\lambda}\left[ \DP{\psi}{\varphi} \left(\Dinv{2}\zeta^2\right)\right]\\ &=& \Dinv{a\cos\varphi}\DP{}{\varphi} \left[\left(\cos\varphi\Dinv{a\cos\varphi}\DP{\psi}{\lambda}\right) \left(\Dinv{2}\zeta^2\right)\right] +\Dinv{a\cos\varphi} \DP{}{\lambda}\left[ -\left(\Dinv{a}\DP{\psi}{\varphi}\right) \left(\Dinv{2}\zeta^2\right)\right]\\ &=& \Dinv{a\cos\varphi}\DP{}{\varphi} \left[\cos\varphi u_\varphi \left(\Dinv{2}\zeta^2\right)\right] +\Dinv{a\cos\varphi} \DP{}{\lambda}\left[ u_\lambda\left(\Dinv{2}\zeta^2\right)\right] \end{eqnarray*} % 第 3 項目は % \begin{eqnarray*} & &\frac{2\Omega}{a^2}\zeta\DP{\psi}{\lambda} = \frac{2\Omega}{a^2} (\Ddiv\Dgrad\psi)\DP{\psi}{\lambda} = \frac{2\Omega}{a^2} \left[ \Ddiv\left(\Dgrad\psi\DP{\psi}{\lambda}\right) -\Dgrad\psi\cdot\DP{\Dgrad\psi}{\lambda} \right]\\ &=& \frac{2\Omega}{a^2} \left[ \Ddiv\left(\Dgrad\psi\DP{\psi}{\lambda}\right) -\frac{1}{2}\DP{(|\Dgrad\psi|^2)}{\lambda}\right]\\ &=& \frac{2\Omega}{a^2} \left\{ \Dinv{a^2\cos\varphi}\DP{}{\lambda}\left( \Dinv{\cos\varphi}\DP{\psi}{\lambda}\DP{\psi}{\lambda}\right) +\Dinv{a^2\cos\varphi}\DP{}{\varphi}\left[\cos\varphi\left( \DP{\psi}{\varphi}\DP{\psi}{\lambda}\right)\right]\right.\\ & & \qquad \left. -\DP{}{\lambda}\frac{1}{2}\left[ \left(\Dinv{a\cos\varphi}\DP{\psi}{\lambda}\right)^2 +\left(\Dinv{a}\DP{\psi}{\varphi}\right)^2\right] \right\}\\ &=& \frac{2\Omega}{a^2} \left\{ \Dinv{a^2\cos\varphi}\DP{}{\lambda} \frac{1}{2}\cos\varphi\left[ \left(\Dinv{\cos\varphi}\DP{\psi}{\lambda}\right)^2 -\left(\DP{\psi}{\varphi}\right)^2 \right] +\Dinv{a^2\cos\varphi}\DP{}{\varphi}\cos\varphi\left( \DP{\psi}{\varphi}\DP{\psi}{\lambda}\right) \right\}\\ &=& \Dinv{a\cos\varphi}\DP{}{\lambda}\left\{ \frac{\Omega}{a}\cos\varphi\left[ \left(\Dinv{a\cos\varphi}\DP{\psi}{\lambda}\right)^2 -\left(\Dinv{a}\DP{\psi}{\varphi}\right)^2 \right]\right\}\\ & & +\Dinv{a^2\cos\varphi}\DP{}{\varphi}\cos\varphi \left(\frac{2\Omega}{a}\cos\varphi \Dinv{a}\DP{\psi}{\varphi} \Dinv{a\cos\varphi}\DP{\psi}{\lambda} \right)\\ &=& \Dinv{a\cos\varphi}\DP{}{\lambda}\left[ \frac{\Omega}{a}\cos\varphi(u_\lambda^2-u_\varphi^2) \right] +\Dinv{a\cos\varphi}\DP{}{\varphi}\left[\cos\varphi\left( \frac{2\Omega}{a}u_\lambda u_\varphi\cos\varphi\right) \right]. \end{eqnarray*} % したがって \begin{eqnarray} \nonumber \DP{}{t}\left(\Dinv{2}\zeta^2\right) &+& \Dinv{a\cos\varphi}\DP{}{\lambda}\left[ u_\lambda\left(\Dinv{2}\zeta^2\right) +\frac{\Omega}{a}\cos\varphi(u_\lambda^2-u_\varphi^2) \right]\\ &+&\Dinv{a\cos\varphi}\DP{}{\varphi}\left\{\cos\varphi \left[ u_\varphi \left(\Dinv{2}\zeta^2\right) + \frac{2\Omega}{a}u_\lambda u_\varphi\cos\varphi\right] \right\} = 0. \Deqlab{相対エンストロフィー保存則} \end{eqnarray} % 相対エンストロフィーの時間変化はフラックスの収束発散という形で 表されている. この式を全球積分すると, 左辺第 2 項経度方向のフラックスの収束は周期境界条件のために 0 となる. 左辺第 3 項緯度方向のフラックスの収束は境界にて $u_\varphi=0$ が 0 であることからこの積分も 0 になる. したがって \begin{equation} \Ddint_S\DP{}{t}\left(\Dinv{2}\zeta^2\right)\Dd S =\DD{}{t}\Ddint_S \Dinv{2}\zeta^2\Dd S =0. \end{equation} % すなわち全球積分した相対エンストフィーは保存する. \subsection{擾乱方程式} \Deqref{非粘性渦度方程式}にあらわれる各変数を 安定性を調べる軸対称な基本場とそれからのずれである擾乱部分に わけて表す. すなわち, % \begin{eqnarray} \Deqlab{振幅展開渦度} \zeta(\lambda,\varphi,t) &=& Z(\varphi) + \zeta'(\lambda,\varphi,t) + \zeta^{(2)} + \cdots,\\ \Deqlab{振幅展開流線} \psi(\lambda,\varphi,t) &=& \Psi(\varphi) + \psi'(\lambda,\varphi,t) + \psi^{(2)} + \cdots,\\ \Deqlab{振幅展開速度経度成分} u_\lambda(\lambda,\varphi,t) &=& U_\lambda(\varphi) + u'_\lambda(\lambda,\varphi,t) + u_\lambda^{(2)} + \cdots,\\ \Deqlab{振幅展開速度緯度成分} u_\varphi(\lambda,\varphi,t) &=& 0 + u'_\varphi(\lambda,\varphi,t) + u_\varphi^{(2)} + \cdots,\\ \end{eqnarray} % ここで $ Z(\varphi), \Psi(\varphi), U_\lambda(\varphi)$ は 安定性を調べる基本場の渦度, 流線, 速度の経度成分である. 基本場は帯状流であるからこれらは経度方向に一様であり, $\varphi$ にのみ依存する. 基本場の速度の緯度成分が 0 であることにも注意されたい. $\zeta',\psi',u'_\lambda, u'_\varphi$ が擾乱の渦度, 流線関数, 速度の経度成分および緯度成分である. 擾乱の大きさは基本場にくらべて十分に小さいものと仮定する. $\zeta^{(2)},\psi^{(2}),u^{(2)}_\lambda, u^{(2)}_\varphi$ は 擾乱 $\zeta',\psi',u'_\lambda, u'_\varphi$ の非線形効果に よって引き起こされるさらに微小な量である. $\zeta',\psi',u'_\lambda, u'_\varphi$ の振幅に対して 大きさがその 2 乗に比例して小さくなるので, $\zeta',\psi',u'_\lambda, u'_\varphi$ を 1 次の微小擾乱, $\zeta^{(2)},\psi^{(2}),u^{(2)}_\lambda, u^{(2)}_\varphi$ を 2 次の微小擾乱とよぶ. \Deqref{振幅展開渦度}, \Deqref{振幅展開流線}を 渦度方程式\Deqref{非粘性渦度方程式}に代入し, 各振幅のオーダーでまとめる. 振幅の 0 次の方程式は左辺が 0 になり, 単に基本場の量が渦度方程式を満たしていることを確認するだけである. 振幅の 1 次の式からは, % \begin{displaymath} \DP{\zeta'}{t} - \Dinv{a^2\cos\varphi}\DP{\Psi}{\varphi}\DP{\zeta'}{\lambda} + \Dinv{a^2\cos\varphi}\DP{\psi}{\lambda}\DP{Z}{\varphi} + \frac{2\Omega}{a^2}\DP{\psi'}{\lambda} = 0, \end{displaymath} % 左辺 2 項目と 3 項目を整理して, 速度に書き直すと \begin{displaymath} \DP{\zeta'}{t} + U_\lambda(\varphi) \Dinv{a\cos\varphi}\DP{\zeta'}{\lambda} + \hat{\beta}(\varphi) \Dinv{a\cos\varphi}\DP{\psi'}{\lambda} = 0, \end{displaymath} % ただし % \begin{equation} \hat{\beta}(\varphi) = \frac{2\Omega\cos\varphi}{a}+ \Dinv{a}\DP{Z}{\varphi} \end{equation} は 基本場のポテンシャル渦度 \begin{equation} Q(\varphi) = 2\Omega\sin\varphi + Z(\varphi) \end{equation} % の緯度微分である. すなわち \begin{equation} \hat{\beta}(\varphi) = \Ddsty \Dinv{a}\DD{Q}{\varphi}. \end{equation} 擾乱の渦度方程式に $\zeta'$ をかけて変形すると, % \begin{displaymath} \DP{}{t}\left(\Dinv{2}\zeta'^2\right) + U_\lambda \Dinv{a\cos\varphi}\DP{}{\lambda} \DP{}{\lambda}\left(\Dinv{2}\zeta'^2\right) + \hat{\beta}(\varphi) \Dinv{a\cos\varphi}\zeta'\DP{\psi'}{\lambda} = 0, \end{displaymath} % $\hat{\beta}(\varphi)/(a\cos\varphi)$ で各辺を割ると % \begin{displaymath} \DP{}{t}\left(\frac{a\cos\varphi}{2\hat{\beta}(\varphi)}\zeta'^2\right) + \Dinv{a\cos\varphi}\DP{}{\lambda}\left[ U_\lambda(\varphi) \left(\frac{a\cos\varphi}{2\hat{\beta}(\varphi)}\zeta'^2\right) \right] + \zeta'\DP{\psi'}{\lambda} = 0, \end{displaymath} % 第 3 項目は\Deqref{相対エンストロフィー保存則}の導出で行なった 式変形を同様に行なうことにより, % \begin{eqnarray*} \zeta'\DP{\psi'}{\lambda} &=& \Dinv{a\cos\varphi}\DP{}{\lambda}\left\{ \left[\Dinv{2} \left(\Dinv{a\cos\varphi}\DP{\psi'}{\lambda}\right)^2 -\left(\Dinv{a}\DP{\psi'}{\varphi}\right)^2 \right]\right\}\\ & & +\Dinv{a^2\cos\varphi}\DP{}{\varphi}\cos\varphi \left(\Dinv{a}\DP{\psi'}{\varphi} \Dinv{a\cos\varphi}\DP{\psi'}{\lambda} \right)\\ &=& \Dinv{a\cos\varphi}\DP{}{\lambda}\left[ \Dinv{2}\cos\varphi(u_\lambda'^2-u_\varphi'^2) \right] +\Dinv{a\cos\varphi}\DP{}{\varphi}\left[\cos\varphi \left(u_\lambda' u_\varphi'\cos\varphi\right) \right]. \end{eqnarray*} % したがって, \begin{eqnarray*} \DP{}{t}\left(\frac{a\cos\varphi}{2\hat{\beta}(\varphi)}\zeta'^2\right) &+& \Dinv{a\cos\varphi}\DP{}{\lambda}\left[ U_\lambda(\varphi) \left(\frac{a\cos\varphi}{2\hat{\beta}(\varphi)}\zeta'^2\right) + \Dinv{2}\cos\varphi(u_\lambda'^2-u_\varphi'^2) \right]\\ &+&\Dinv{a\cos\varphi}\DP{}{\varphi}\left[\cos\varphi \left(u_\lambda' u_\varphi'\cos\varphi\right) \right] = 0, \end{eqnarray*} % ここであらためて \begin{displaymath} {\cal A} \equiv \frac{a\cos\varphi}{2\hat{\beta}(\varphi)}\zeta'^2, \quad {\cal F}_\lambda \equiv U_\lambda(\varphi) {\cal A} + \Dinv{2}\cos\varphi(u_\lambda'^2-u_\varphi'^2), \quad {\cal F}_\varphi \equiv u_\lambda' u_\varphi'\cos\varphi, \end{displaymath} % と表せば \begin{equation} \DP{\cal A}{t} + + \Dinv{a\cos\varphi}\DP{{\cal F}_\lambda}{\lambda} + \Dinv{a\cos\varphi}\DP{}{\varphi}({\cal F}_\varphi\cos\varphi) =0. \end{equation} % 1 次擾乱の 2 乗の量(振幅の 2 次のオーダー) $A$ は フラックス形式の保存則にしたがうことがわかる. すなわち, 全球で積分すればフラックス $({\cal F}_\lambda,{\cal F}_\varphi)$ の収束発散項は 0 となり \begin{equation} \Deqlab{1次擾乱の2乗の保存則} \Ddint_S \DP{\cal A}{t} \Dd S = \DD{}{t}\Ddint_S {\cal A} \Dd S =0, \end{equation} % したがって ${\cal A}$ は保存量 \footnote{天下り的に導出した保存量 ${\cal A}$ は, 実は 2 次のオーダーの角運動量 $u^{(2)}a\cos\varphi$ と 関連づけることができる(ここでは示さない). } である. \subsection{不安定の発生の必要条件(Rayleigh's criterion)} 1 次擾乱の 2 乗の保存則\Deqref{1次擾乱の2乗の保存則}から, 不安定が発生するための基本場に関する必要条件を導くことができる. \Deqref{1次擾乱の2乗の保存則}を改めて書くと, % \begin{equation} \DD{}{t}\Ddint_S {\cal A} \Dd S =\DD{}{t}\Ddint_S \frac{a\cos\varphi}{2\hat{\beta}(\varphi)}\zeta'^2 \Dd S =\DD{}{t}\Ddint_S \frac{a^2\cos^2\varphi}{2\hat{\beta}(\varphi)}\zeta'^2 \Dd \lambda \Dd \varphi =0, \end{equation} % である. 基本場が不安定であると擾乱の振幅が増加するので $\zeta'^2$ も時間的に増加していく. そのような状況でも上の保存則は満たされなければならない. そのためには, 保存量 ${\cal A}$ は 0 でなければならない. % \begin{equation} \Ddint_S \frac{a^2\cos^2\varphi}{2\hat{\beta}(\varphi)}\zeta'^2 \Dd \lambda \Dd \varphi =0, \end{equation} % すなわち, 被積分関数の分布が正のところと負のところが存在して, それぞれの振幅が増大していくが, 積分値は常にキャンセルして 0 となっている. 被積分関数の正負は $\zeta'^2$ の係数の $\hat{\beta}(\varphi)$ で 定められることになる. したがって基本場が不安定となるためには, 基本場のポテンシャル渦度の緯度方向の傾度 $\hat{\beta}(\varphi)$ が 符号を変えねばならない. いいかえれば % \begin{equation} \hat{\beta}(\varphi) = \frac{2\Omega\cos\varphi}{a}+ \Dinv{a}\DP{Z}{\varphi} = 0, \end{equation} となる緯度が存在することが必要である. これが基本場が不安定となるのための必要条件である. \section{実験設定} \subsection{初期条件} 安定性を調べる基本流は, その渦度分布が帯状波数 0 の球面調和函数 1 成分で表される帯状流である. % \begin{eqnarray} \zeta_b &=& \zeta_0 P_n^0(\cos\varphi), \\ \psi_b &=& \nabla^{-2}\zeta_b = -\frac{a^2\zeta_0}{n(n+1)} P_n^0(\cos\varphi),\\ u_{\lambda,b} &=& -\frac{1}{a}\DP{\psi_b}{\varphi} = \frac{a\zeta_0}{n(n+1)}\DP{P_n^0(\cos\varphi)}{\varphi},\\ u_{\varphi,b} &=& \frac{1}{a\cos\varphi}\DP{\psi_b}{\lambda} =0, \end{eqnarray} % ここで, 下付き添字 $b$ は基本場の量であることを示す. この基本流に対して微小振幅のランダムな渦度のノイズを加えたものを 初期条件とする. \subsection{パラメター} 問題に登場する物理パラメターは \begin{itemize} \item 球の半径 $a$ \item 帯状流の全波数 $n$ \item 基本場の相対渦度と系の回転角速度の比 $\zeta_0/\Omega$ \item ノイズの振幅 \item ランダムノイズの生成のためのパラメター \end{itemize} % このうち, 長さスケールを球の半径にとることにして $a=1$ と することができる. 上記の他に, 粘性過程を定めるパラメターがある. % \begin{itemize} \item 高階粘性拡散の次数 \item 高階粘性拡散係数 \end{itemize} % 理想的には粘性のない状態での解を知りたいのであるが, 数値計算の安定性のためには高波数成分エネルギーをおさえる 必要があるので, 高階粘性拡散を導入している. したがって以下のパラメターは, 知りたい現象の波数においては時間積分する計算時間にわたって減衰しない 程度にし, 逆に切断波数近辺の高波数成分に対しては計算時間内に 粘性により充分に減衰させるようにチューニングすべきである. % \subsection{時間積分} % 数値計算が意味を持つための時間積分の時間ステップの目安を見積もる. % 各時間変化項の最大値は % % % \begin{eqnarray*} % \Dinv{a^2}J(\psi,\zeta) % & \sim & \frac{UL\cdot U/L}{L^2} % \sim \frac{U^2}{L^2},\\ % \frac{2\Omega}{a^2}\DP{\psi}{\lambda} % & \sim & \frac{\zeta}{a^2}\cdot Ua % \sim \frac{\zeta U}{a},\\ % (-1)^{p+1}\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p \zeta % & \sim & \nu_{2p}\left(\frac{N_{max}(N_{max}+1)}{a^2}\right)^p % \cdot\frac{U}{L} % \end{eqnarray*} % % % ここで, $U,L$ は現象の典型的な速度と水平スケールである. % 一方, 時間変化項は % $\Ddsty \DP{\zeta}{t}\sim \frac{U}{(L)(\Delta t)}$ % と見積もられるので, 時間変化が十分に解像できるには % $\Ddsty \DP{\zeta}{t}$ が各時間変化項の見積もりより大きくなるよう % 時間刻 $\Delta t$ をとらねばならない % \footnote{どうもすっきりしない論理だ...}. % したがって, % % % \begin{eqnarray*} % \frac{U}{L\Delta t} % & \gg & \max\left[ % \Dinv{a^2}J(\psi,\zeta), % \frac{2\Omega}{a^2}\DP{\psi}{\lambda}, % (-1)^{p+1}\nu_{2p} % \left(\Dlapla + \frac{2}{a^2}\right)^p \zeta % \right]\\ % &\sim &\max\left[ % \frac{U^2}{L^2}, % \frac{\zeta U}{a}, % \nu_{2p}\left(\frac{N_{max}(N_{max}+1)}{a^2}\right)^p % \cdot\frac{U}{L} % \right],\\ % \Delta t % &\ll& \min\left[ % \frac{L}{U}, % \frac{a}{\zeta L}, % \Dinv{\nu_{2p}\left(\frac{N_{max}(N_{max}+1)}{a^2}\right)^p} % \right]. % \end{eqnarray*} % % % 右辺の各項が最小となる水平スケールをそれぞれ選べば, % 第 1 項目は $L \sim \Delta x$, 第 2 項目は $L \sim \Delta a$ % となる. したがって, % \begin{equation} % \Delta t % \ll \min\left[ % \frac{\Delta x}{U}, % \frac{1}{\Omega}, % \Dinv{\nu_{2p}\left(\frac{N_{max}(N_{max}+1)}{a^2}\right)^p} % \right]. % \end{equation} \subsection{課題} \begin{enumerate} \item 剛体回転流($n=1$)のみをノイズなしの初期条件としてあたえ, それが振幅によらず定常に保たれることを確認せよ. \item 剛体回転流($n=1$)にノイズを与えた場合に, 回転角速度と基本流の渦度の比をいろいろ変えてみて 基本流が安定であることを確認せよ. \item $n=2$ の基本流についても同様に行なえ. \item $n=3$ 以上の波数の基本流は不安定となる可能性がある. Bains (1976) の表 1 と図 1 を参考にして, $n=3$ から順番に, 不安定となると予想される基本流の渦度と回転角速度との比の値より 大きい値を与えた場合と小さい値を与えた場合の時間発展をそれぞれ 計算し, 臨界値をこえた状況でのみ基本流が不安定となることを確認せよ. また, 発生する不安定場の帯状波数が予想されているものと一致しているか も確認せよ. \item 経度平均したポテンシャル渦度の緯度傾度分布 $\hat{\beta}(\varphi)$ を描き, 不安定が発生する状況では基本場(初期状態)の $\hat{\beta}(\varphi)$ が 0 となる点(広義の変曲点)が 存在していることを確かめよ. また, 主に変曲点周辺にて擾乱が発生していることを確かめよ. \item 各スキームでの数値的に安定に計算するための時間刻について議論せよ. \end{enumerate} \begin{Dreference} \item SPMODEL 解説文書「2 次元回転球面上での非発散順圧流体の定式化」 http://www.gfd-dennou.org/library/spmodel/2d-sphere-w/baro/pub/ \item Baines, P. G., 1976 : The stability of planetary waves on a sphere. {\it J. Fluid Mech.}, {\bf 73}, 193--223. \end{Dreference} \end{document} %%% Local Variables: %%% mode: latex %%% TeX-master: t %%% End: