% 表題 回転球面上での 2 次元順圧流体の自由減衰乱流問題の定式化 % % 履歴 2005/03/09 竹広 真一 % 2007/11/25 竹広 真一 参照文書追加 % 2008/01/03 竹広 真一 エネルギースペクトル計算追加 % 2008/02/12 竹広 真一 ロスビー波振動数バグfix % 2008/02/25 竹広 真一 全運動エネルギー, 全エンストロフィー時間変化追加 % 2008/03/02 竹広 真一 全運動エネルギー, 全エンストロフィー時間変化修正 % \documentclass[12pt,a4j,notitlepage]{jarticle} \usepackage{Dennou6} \usepackage{psfrag} \usepackage{ascmac} \usepackage{times} \Dtitle[回転球面 2 次元順圧自由減衰乱流] {回転球面上での 2 次元順圧流体の\\自由減衰乱流問題の定式化} \Dauthor[竹広 真一]{竹広 真一, SPMODEL 開発グループ} \Ddate{平成 20 年 3 月 4 日} \Dparskip \pagestyle{Dheadings} \begin{document} \maketitle \section{はじめに} \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{定義} 球面上の各点$(\lambda,\varphi)$における局所的な 運動エネルギー $e_k$ は次のように定義される. % \begin{equation} \Deqlab{局所運動エネルギー} e_k(\lambda,\varphi,t)=\frac{1}{2}|\Dvect{u}|^2 =\frac{1}{2}(u_\lambda^2+u_\varphi^2). \end{equation} % 全運動エネルギー $E_k$ はこれを球面全体で積分して得られる. \begin{eqnarray*} \Deqlab{全運動エネルギー} E_k(t) &=& \Ddint_S\Dd S e_k(\lambda,\varphi,t) = \Dinv{2}\Ddint_S\Dd S |\Dvect{u}|^2 \\ &=& \Dinv{2}\int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \ a^2\cos\varphi (u_\lambda^2+u_\varphi^2). \end{eqnarray*} 球面上の各点$(\lambda,\varphi)$における局所的な エンストロフィーは $\Ddsty \frac{1}{2}\zeta^2$ で定義される. 全エンストロフィー $\cal Q$ はこれを球面全体で積分して得られる. % \begin{equation} \Deqlab{全エンストロフィー} {\cal Q}(t) = \Ddint_S \Dd S \left(\frac{1}{2}\zeta^2\right) = \Dinv{2}\int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \ a^2\cos\varphi\ \zeta^2 \end{equation} \subsection{エネルギー・エンストロフィーのスペクトル分解} \Deqref{速度流線関数関係}を用いて全エネルギーを表現すると % \begin{eqnarray*} E_k(t) &=& \Dinv{2}\int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \ a^2\cos\varphi \left[ \left(\Dinv{a}\DP{\psi}{\varphi}\right)^2 + \left(\Dinv{a\cos\varphi}\DP{\psi}{\lambda}\right)^2 \right] \\ &=& \Dinv{2}\int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \ \cos\varphi \left[ \left(\DP{\psi}{\varphi}\right)^2 + \left(\Dinv{\cos\varphi}\DP{\psi}{\lambda}\right)^2 \right]. \end{eqnarray*} % ここで流線関数を球面調和函数で展開する. すなわち, % \begin{equation} \psi(\lambda,\varphi) =\sum_{n=0}^\infty\sum_{m=-n}^n \tilde{\psi}_n^m(t)Y_n^m(\lambda,\varphi) \end{equation} % これを上の式に代入すると % \begin{eqnarray*} E_k(t) &=& \Dinv{2}\int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \ \cos\varphi \left[ \left(\sum_{n=0}^\infty\sum_{m=-n}^n \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \tilde{\psi}_n^m\DP{Y_n^m}{\varphi} \tilde{\psi}_{n'}^{m'} \DP{Y_{n'}^{m'}}{\varphi}\right) \right.\\ & &\left.\hspace{35mm} + \left(\Dinv{\cos^2\varphi} \sum_{n=0}^\infty\sum_{m=-n}^n \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \tilde{\psi}_n^m\DP{Y_n^m}{\lambda} \tilde{\psi}_{n'}^{m'} \DP{Y_{n'}^{m'}}{\lambda}\right) \right]\\ &=& \Dinv{2}\sum_{n=0}^\infty\sum_{m=-n}^n \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'}\\ & & \hspace*{25mm} \int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \ \cos\varphi \left( \DP{Y_n^m}{\varphi}\DP{Y_{n'}^{m'}}{\varphi} +\Dinv{\cos^2\varphi} \DP{Y_n^m}{\lambda}\DP{Y_{n'}^{m'}}{\lambda} \right)\\ &=& \Dinv{2}\sum_{n=0}^\infty\sum_{m=-n}^n \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'} \int_0^{2\pi}\Dd \lambda \left[Y_n^m\DP{Y_{n'}^{m'}}{\varphi}\cos\varphi \right]_{-\pi/2}^{\pi/2}\\ & & \hspace*{45mm} -\int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \ Y_n^m\DP{}{\varphi}\cos\varphi \DP{Y_{n'}^{m'}}{\varphi}\\ & & \hspace*{45mm} +\int_{-\pi/2}^{\pi/2}\Dd \varphi \ \left[ \Dinv{\cos^2\varphi} Y_n^m\DP{Y_{n'}^{m'}}{\lambda} \right]_0^{2\pi}\\ & & \hspace*{45mm} -\int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \ \Dinv{\cos^2\varphi} Y_n^m\DP[2]{Y_{n'}^{m'}}{\lambda} \\ &=& -\Dinv{2}\sum_{n=0}^\infty\sum_{m=-n}^n \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'}\\ & & \hspace*{10mm} \int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \cos\varphi Y_n^m \left[\Dinv{\cos\varphi} \DP{}{\varphi} \left( \cos\varphi \DP{Y_{n'}^{m'}}{\varphi} \right) + \Dinv{\cos^2\varphi} \DP[2]{Y_{n'}^{m'}}{\lambda} \right]\\ &=& \Dinv{2}\sum_{n=0}^\infty\sum_{m=-n}^n \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'} \int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \cos\varphi \ n'(n'+1)Y_n^m Y_{n'}^{m'}\\ &=& \Dinv{2}\sum_{n=0}^\infty\sum_{m=-n}^n \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'} \delta_{nn'}\delta_{mm'} n'(n'+1) N_l^m\\ &=& \Dinv{2}\sum_{n=0}^\infty\sum_{m=-n}^n N_l^m n(n+1) |\tilde{\psi}_n^m(t)|^2. \end{eqnarray*} % ここで $N_n^m$ は % \begin{equation} N_n^m \equiv \Ddint_S \Dd S |Y_n^m(\lambda,\varphi)|^2 = \int_0^{2\pi}\Dd\lambda\int_{-\pi/2}^{\pi/2} \Dd\varphi \cos\varphi|Y_n^m(\lambda,\varphi)|^2, \end{equation} % であり, 球面調和函数の規格化によって値が異なる \footnote{ISPACK では 2 に規格化されたルジャンドル函数を用いているので $N_n^m=4\pi$ である(要確認)}. したがって, エネルギースペクトル密度 ${\mathcal E}(n,t)$を % \begin{equation} \Deqlab{エネルギースペクトル密度} {\mathcal E}(n,t) \equiv \Dinv{2}\sum_{m=-n}^n N_n^m n(n+1) |\tilde{\psi}_n^m(t)|^2 \end{equation} % と定義できて, % \begin{equation} E_k(t)=\sum_{n=0}^\infty {\mathcal E}(n,t), \end{equation} % となる. 特に, 球面調和函数が $4\pi$ で規格化されている場合, すなわち $N_n^m=4\pi$ の場合には, この規格化係数を省略したエネルギー密度を % \begin{equation} \Deqlab{エネルギースペクトル密度係数省略} {\mathcal E}(n,t) \equiv \Dinv{2}n(n+1)\sum_{m=-n}^n |\tilde{\psi}_n^m(t)|^2 \end{equation} % と定義することができる. その場合, % \begin{equation} E_k(t)=4 \pi \sum_{n=0}^\infty {\mathcal E}(n,t), \end{equation} % となる. 同様に, \Deqref{渦度流線関数関係}を用いて全エンストロフィーを表し, 流線関数を球面調和函数展開すれば, % \begin{eqnarray*} {\cal Q}(t) &=& \frac{1}{2}\int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \ a^2\cos\varphi\ (\Dlapla \psi)^2\\ &=& \frac{1}{2}\int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \ a^2\cos\varphi\ \left(\Dlapla\sum_{n=0}^\infty\sum_{m=-n}^n \tilde{\psi}_n^m Y_n^m \right) \left(\Dlapla \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \tilde{\psi}_{n'}^{m'}Y_{n'}^{m'} \right)\\ &=& \frac{1}{2} \sum_{n=0}^\infty\sum_{m=-n}^n \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \ a^2\cos\varphi\ \\ & & \hspace{50mm} \times\left( -\frac{n(n+1)}{a^2}\tilde{\psi}_n^m Y_n^m \right) \left(-\frac{n'(n'+1)}{a^2} \tilde{\psi}_{n'}^{m'}Y_{n'}^{m'} \right)\\ &=& \frac{1}{2} \sum_{n=0}^\infty\sum_{m=-n}^n \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \frac{n(n+1)n'(n'+1)}{a^2} \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'} \int_0^{2\pi}\Dd \lambda \int_{-\pi/2}^{\pi/2}\Dd \varphi \ \cos\varphi\ Y_n^m Y_{n'}^{m'}\\ &=& \frac{1}{2}\sum_{n=0}^\infty\sum_{m=-n}^n \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \frac{n(n+1)n'(n'+1)}{a^2} \tilde{\psi}_n^m\tilde{\psi}_{n'}^{m'} N_l^m \delta_{n,n'} \delta_{m,m'}\\ &=& \frac{1}{2}\sum_{n=0}^\infty\sum_{m=-n}^n N_l^m \frac{n^2(n+1)^2}{a^2}|\tilde{\psi}_n^m(t)|^2. \end{eqnarray*} % したがって, エンストロフィースペクトル密度 ${\mathcal Q}(n,t)$を % \begin{equation} \Deqlab{エンストロフィースペクトル密度} {\mathcal Q}(n,t) \equiv \Dinv{2a^2}\sum_{m=-n}^n N_n^m n^2(n+1)^2 |\tilde{\psi}_n^m(t)|^2 \end{equation} % と定義できて, % \begin{equation} Q(t)=\sum_{n=0}^\infty {\mathcal Q}(n,t), \end{equation} % となる. エンストロフィーについても, 球面調和函数が $4\pi$ で規格化されている場合, すなわち $N_l^m=4\pi$ の場合には, 規格化係数を省略したエンストロフィー密度を % \begin{equation} \Deqlab{エンストロフィースペクトル密度係数省略} {\mathcal Q}(n,t) \equiv \frac{n^2(n+1)^2}{2a^2}\sum_{m=-n}^n |\tilde{\psi}_n^m(t)|^2 \end{equation} % と定義することができる. その場合, % \begin{equation} Q(t)=4 \pi \sum_{n=0}^\infty {\mathcal Q}(n,t), \end{equation} % となる. \subsection{全運動エネルギー, 全エンストロフィーの時間変化} 全運動エネルギーの式は渦度方程式\Deqref{渦度方程式}に $-\psi$ を かけて球面全体で積分することにより得ることができる. % \begin{displaymath} -\int\int dS \psi\DP{\zeta}{t} -\int\int dS \psi\Dinv{a^2}J(\psi,\zeta) -\int\int dS \psi\frac{2\Omega}{a^2}\DP{\psi}{\lambda} = -\int\int dS \psi(-1)^{p+1}\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p \zeta, \end{displaymath} % ただし % \begin{displaymath} \int\int dS =a^2\int_{0}^{2\pi}d\lambda\ \int_{-\pi/2}^{\pi/2}d\varphi\cos\varphi = a^2\int_{0}^{2\pi}d\lambda\int_{-1}^{1}d\mu \end{displaymath} % は球面全積分である. 第 1 項目は \begin{eqnarray*} & & -\int\int dS \psi\DP{\zeta}{t} = -\int\int dS \psi\DP{}{t} \Dlapla \psi\\ & & = -\int\int dS \psi\DP{}{t} \left[ \Dinv{a^2\cos^2\varphi}\DP[2]{\psi}{\lambda} + \Dinv{a^2\cos\varphi}\DP{}{\varphi} \left(\cos\varphi\DP{\psi}{\varphi}\right)\right]\\ & & = \int\int dS \Dinv{a^2\cos^2\varphi}\DP{\psi}{\lambda} \Dinv{a^2\cos^2\varphi}\DP{}{\lambda} \DP{\psi}{t} + \int\int dS \Dinv{a^2}\DP{\psi}{\varphi}\DP{}{\varphi}\DP{\psi}{t}\\ & & = \int\int dS \frac{1}{2}\DP{}{t}\Dinv{a^2\cos^2\varphi}\left(\DP{\psi}{\lambda}\right)^2 + \int\int dS \frac{1}{2}\DP{}{t} \Dinv{a^2}\left(\DP{\psi}{\varphi}\right)^2\\ & & = \DD{}{t}\int\int dS \frac{1}{2}(u_\lambda^2+u_\varphi^2) = \DD{E_k}{t} \end{eqnarray*} % 第 2 項目は \begin{eqnarray*} \int\int dS \psi\Dinv{a^2}J(\psi,\zeta) = \int\int dS \Dinv{a^2}J(\frac{1}{2}\psi^2,\zeta) \end{eqnarray*} % ここでヤコビアン $J(f,g)$ は % \begin{eqnarray*} J(f,g) = \DP{f}{\lambda}\DP{g}{\mu}-\DP{f}{\mu}\DP{g}{\lambda} = \DP{}{\lambda}\left(f\DP{g}{\mu}\right) -\DP{}{\mu}\left(f\DP{g}{\lambda}\right) \end{eqnarray*} % と変形できるので全球積分すると 0 となる. % 第 3 項目は \begin{eqnarray*} -\int\int dS \psi\frac{2\Omega}{a^2}\DP{\psi}{\lambda} = -\int\int dS \frac{\Omega}{a^2}\DP{\psi^2}{\lambda} = 0. \end{eqnarray*} % 最後の粘性項は球面調和函数展開して評価する. % \begin{eqnarray*} & & -\int\int dS \psi(-1)^{p+1}\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p \zeta = \int\int dS \psi (-1)^p\nu_{2p}\left(\Dlapla + \frac{2}{a^2}\right)^p \Dlapla\psi\\ & & = \int\int dS \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \tilde{\psi}_{n'}^{m'}(t)Y_{n'}^{m'}(\lambda,\varphi)\\ & &\qquad\qquad \sum_{n=0}^\infty\sum_{m=-n}^n (-1)^p\nu_{2p} \left(\frac{-n(n+1) + 2}{a^2}\right)^p \left(\frac{-n(n+1)}{a^2}\right) \tilde{\psi}_n^m(t)Y_n^m(\lambda,\varphi) \\ & & = -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p \left(\frac{n(n+1)}{a^2}\right) |\tilde{\psi}_n^m(t)|^2 \int\int dS Y_n^m(\lambda,\varphi)^2 \\ & & = -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p \left(\frac{n(n+1)}{a^2}\right) a^2 N_n^m |\tilde{\psi}_n^m(t)|^2\\ & & = -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p n(n+1) N_n^m |\tilde{\psi}_n^m(t)|^2 \end{eqnarray*} % ただし, $N_n^m$ は球面調和函数の規格化定数である. よって全運動エネルギーの時間変化は % \begin{equation} \DD{E_k}{t} = -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p n(n+1) N_n^m |\tilde{\psi}_n^m(t)|^2. \end{equation} % 特に球面調和函数の規格化定数が $N_n^m=4\pi$ の場合には % \begin{eqnarray} \DD{E_k}{t} &=& -4\pi\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p n(n+1) |\tilde{\psi}_n^m(t)|^2 \\ &=& -4\pi\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p 2 {\mathcal E}(n,t). \end{eqnarray} % ただし ${\mathcal E}(n,t)$ は\Deqref{エネルギースペクトル密度係数省略}で定義さ れるスペクトル密度である. 全運動エネルギーの代わりに球面平均したエネルギー $\bar{E_k}$ の時間変化は 全体を $4\pi a^2$ でわって % \begin{eqnarray} \DD{\bar{E_k}}{t} &=& -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p \frac{n(n+1)}{a^2} |\tilde{\psi}_n^m(t)|^2\\ &=& -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p \frac{2}{a^2} {\mathcal E}(n,t). \end{eqnarray} % $E_k, \bar{E_k}$ も球面調和函数展開で表しておくと % \begin{eqnarray} E_k &=& 4\pi\sum_{n=0}^\infty\sum_{m=-n}^n \frac{n(n+1)}{2} |\tilde{\psi}_n^m(t)|^2 = 4\pi\sum_{n=0}^\infty {\mathcal E}(n,t), \\ \bar{E_k} &=& \sum_{n=0}^\infty\sum_{m=-n}^n \frac{n(n+1)}{2a^2} |\tilde{\psi}_n^m(t)|^2 = \Dinv{a^2}\sum_{n=0}^\infty {\mathcal E}(n,t). \end{eqnarray} 全エンストロフィーの時間変化の式は渦度方程式\Deqref{渦度方程式}に $\Ddsty \zeta=\Dinv{a^2}\Dlapla\psi$ をかけて 球面全体で積分することにより得ることができる. % \begin{displaymath} \int\int dS \zeta\DP{\zeta}{t} +\int\int dS \zeta\Dinv{a^2}J(\psi,\zeta) \int\int dS \zeta\frac{2\Omega}{a^2}\DP{\psi}{\lambda} = \int\int dS \zeta(-1)^{p+1}\nu_{2p} \left(\Dlapla + \frac{2}{a^2}\right)^p \zeta, \end{displaymath} % 第 1 項目は % \begin{displaymath} \int\int dS \zeta\DP{\zeta}{t} = \int\int dS \DP{}{t}\left(\frac{1}{2}\zeta^2\right) = \DD{Q}{t}. \end{displaymath} % 第 2 項目は \begin{displaymath} \int\int dS \zeta\Dinv{a^2}J(\psi,\zeta) \int\int dS \Dinv{2a^2}J(\psi,\zeta^2) = 0. \end{displaymath} % 第 3 項目は % \begin{eqnarray*} \int\int dS \zeta\frac{2\Omega}{a^2}\DP{\psi}{\lambda} = \frac{2\Omega}{a^2}\int\int dS \left[ \Dinv{a^2\cos^2\varphi}\DP[2]{\psi}{\lambda} + \Dinv{a^2\cos\varphi}\DP{}{\varphi} \left(\cos\varphi\DP{\psi}{\varphi}\right) \right] \DP{\psi}{\lambda} \end{eqnarray*} % 1 つ目の被積分関数は % \begin{eqnarray*} \Dinv{a^2\cos^2\varphi}\DP[2]{\psi}{\lambda} \DP{\psi}{\lambda} = \Dinv{2a^2\cos^2\varphi}\DP[]{}{\lambda} \DP[][]{\psi}{\lambda}^2 \end{eqnarray*} % と変形できて, この全球積分は 0 となる. 2 つ目の被積分関数は % \begin{eqnarray*} & & \Dinv{a^2\cos\varphi}\DP{}{\varphi} \left(\cos\varphi\DP{\psi}{\varphi}\right)\DP{\psi}{\lambda} \\ & & = \DP{}{\lambda}\left[\psi \Dinv{a^2\cos\varphi}\DP{}{\varphi} \left(\cos\varphi\DP{\psi}{\varphi}\right)\right] - \psi\Dinv{a^2\cos\varphi}\DP{}{\varphi} \left(\cos\varphi\frac{\partial^2 \psi} {\partial\varphi\partial\lambda}\right), \end{eqnarray*} % 一方で \begin{eqnarray*} & & \Dinv{a^2\cos\varphi}\DP{}{\varphi} \left(\cos\varphi\DP{\psi}{\varphi}\right)\DP{\psi}{\lambda} = \Dinv{a^2\cos\varphi}\DP{}{\varphi}\left[ \DP{\psi}{\lambda} \left(\cos\varphi\DP{\psi}{\varphi}\right)\right] - \frac{\cos\varphi}{a^2\cos\varphi}\DP{\psi}{\varphi} \frac{\partial^2 \psi} {\partial\varphi\partial\lambda}\\ & & = \Dinv{a^2\cos\varphi}\DP{}{\varphi}\left[ \DP{\psi}{\lambda} \left(\cos\varphi\DP{\psi}{\varphi}\right)\right] - \Dinv{a^2\cos\varphi}\DP{}{\varphi}\left[ \cos\varphi\psi \frac{\partial^2 \psi} {\partial\varphi\partial\lambda}\right]\\ & & \qquad + \psi\Dinv{a^2\cos\varphi}\DP{}{\varphi}\left[ \cos\varphi \frac{\partial^2 \psi} {\partial\varphi\partial\lambda}\right]. \end{eqnarray*} % これらを1/2 倍して足し合わせると最後の項がキャンセルして % \begin{eqnarray*} & & \Dinv{a^2\cos\varphi}\DP{}{\varphi} \left(\cos\varphi\DP{\psi}{\varphi}\right)\DP{\psi}{\lambda}\\ & & = \Dinv{2}\DP{}{\lambda}\left[\psi \Dinv{a^2\cos\varphi}\DP{}{\varphi} \left(\cos\varphi\DP{\psi}{\varphi}\right)\right] + \Dinv{2a^2\cos\varphi}\DP{}{\varphi}\left[ \DP{\psi}{\lambda} \left(\cos\varphi\DP{\psi}{\varphi}\right)\right]\\ & & \qquad - \Dinv{2a^2\cos\varphi}\DP{}{\varphi}\left[ \cos\varphi\psi \frac{\partial^2 \psi} {\partial\varphi\partial\lambda}\right]. \end{eqnarray*} % これを全球積分すると 0 になる. 最後に粘性項を球面調和関数展開で表現すると % \begin{eqnarray*} & & \int\int dS \zeta (-1)^{p+1}\nu_{2p} \left(\Dlapla + \frac{2}{a^2}\right)^p \zeta = \int\int dS \Dlapla\psi (-1)^{p+1}\nu_{2p} \left(\Dlapla + \frac{2}{a^2}\right)^p \Dlapla\psi\\ & & = \int\int dS (-1)^{p+1}\nu_{2p} \sum_{n'=0}^\infty\sum_{m'=-n'}^{n'} \left(\frac{-n'(n'+1)}{a^2}\right) \tilde{\psi}_{n'}^{m'}(t)Y_{n'}^{m'}(\lambda,\varphi)\\ & & \qquad \sum_{n=0}^\infty\sum_{m=-n}^n (-1)^{p+1}\nu_{2p} \left(\frac{-n(n+1) + 2}{a^2}\right)^p \left(\frac{-n(n+1)}{a^2}\right) \tilde{\psi}_n^m(t)Y_n^m(\lambda,\varphi) \\ & & = \sum_{n=0}^\infty\sum_{m=-n}^n (-1)^{p+1}\nu_{2p} \left(\frac{-n(n+1) + 2}{a^2}\right)^p \left(\frac{-n(n+1)}{a^2}\right)^2 |\tilde{\psi}_n^m(t)|^2 \int\int dS |Y_n^m(\lambda,\varphi)|^2 \\ & & = -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p \frac{n^2(n+1)^2}{a^2}N_n^m|\tilde{\psi}_n^m(t)|^2 \\ & & = -\sum_{n=0}^\infty \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p 2 {\mathcal Q}(n,t) \end{eqnarray*} % したがって全エンストロフィーの時間変化の式は % \begin{equation} \DD{Q}{t} = -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p \frac{n^2(n+1)^2}{a^2}N_n^m|\tilde{\psi}_n^m(t)|^2\\ \end{equation} % 特に球面調和函数の規格化係数 $N_n^m=4\pi$ の場合には 全エンストロフィーの時間変化の式は % \begin{eqnarray*} \DD{Q}{t} &=& -4\pi\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p \frac{n^2(n+1)^2}{a^2}|\tilde{\psi}_n^m(t)|^2\\ &=& -4\pi\sum_{n=0}^\infty \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p 2 {\mathcal Q}(n,t). \end{eqnarray*} % ただし ${\mathcal Q}(n,t)$ は\Deqref{エンストロフィースペクトル密度係数省略}で 定義されるエンストロフィースペクトル密度である. 全エンストロフィーの代わりに球面平均したエンストロフィー $\bar{Q}$ の 時間変化は全体を $4\pi a^2$ でわって % \begin{eqnarray*} \DD{\bar{Q}}{t} &=& -\sum_{n=0}^\infty\sum_{m=-n}^n \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p \frac{n^2(n+1)^2}{a^4}|\tilde{\psi}_n^m(t)|^2\\ &=& -\sum_{n=0}^\infty \nu_{2p} \left(\frac{n(n+1) - 2}{a^2}\right)^p \frac{2}{a^2} {\mathcal Q}(n,t). \end{eqnarray*} % $Q,\bar{Q}$ も球面調和関数展開表現で表しておくと \begin{eqnarray} Q &=& 4\pi\frac{1}{2}\sum_{n=0}^\infty\sum_{m=-n}^n \frac{n^2(n+1)^2}{a^2}|\tilde{\psi}_n^m(t)|^2 = 4\pi\sum_{n=0}^\infty {\mathcal Q}(n,t), \\ \bar{Q} &=& \frac{1}{2}\sum_{n=0}^\infty\sum_{m=-n}^n \frac{n^2(n+1)^2}{a^4}|\tilde{\psi}_n^m(t)|^2 = \frac{1}{a^2}\sum_{n=0}^\infty {\mathcal Q}(n,t), \\ \end{eqnarray} \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[\frac{-n(n+1)+ 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[\frac{-n(n+1)+2}{a^2}\right]^p = \nu_{2p}\left[\frac{n(n+1)-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 = -a^2\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}{n(n+1)} \tilde{\zeta}_n^{-m},\quad \frac{2\Omega}{a^2}\left[\DP{\psi}{\lambda}\right]_n^{-m} = -\frac{2\Omega m}{ 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}{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}{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|}{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{初期条件} \subsection{パラメター} \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{\Omega}{a^2}\cdot Ua \sim \frac{\Omega 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{\Omega 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}{\Omega 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 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} \begin{Dreference} \item SPMODEL 解説文書「2 次元回転球面上での非発散順圧流体の定式化」 http://www.gfd-dennou.org/library/spmodel/2d-sphere-w/baro/pub/ \end{Dreference} \end{document} %%% Local Variables: %%% mode: latex %%% TeX-master: t %%% End: