\underline{回転球殻 MHD ダイナモの定式化目次 1}

回転球殻 MHD ダイナモの定式化

佐々木洋平, 竹広真一, 林祥介

平成 14 年 10 月 4 日

目次

1 基礎方程式系 2 1.1 系の設定と支配方程式 . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 ポテンシャルによる方程式系の表現 . . . . . . . . . . . . . . . . . . 6

1.2.1 速度場及び磁場とポテンシャルとの関係式 . . . . . . . . . . 7 1.3 境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

1.3.1 運動学的境界条件 . . . . . . . . . . . . . . . . . . . . . . . 9 1.3.2 力学的境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.3.3 熱的境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.3.4 磁場の境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . 10

2 水平離散化 11 2.1 スペクトル展開 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2 方程式のスペクトル表現 . . . . . . . . . . . . . . . . . . . . . . . . 13

2.2.1 運動方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.2 熱の式 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.3 誘導方程式 . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.2.4 速度場の表現 . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.5 磁場の表現 . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3 鉛直離散化 15 3.1 鉛直離散化 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 3.2 鉛直境界条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

3.2.1 トロイダルポテンシャル . . . . . . . . . . . . . . . . . . . . 16 3.2.2 ポロイダルポテンシャル . . . . . . . . . . . . . . . . . . . . 16 3.2.3 温度 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.2.4 磁場 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 \underline{回転球殻 MHD ダイナモの定式化目次 2} 参考文献 20

\underline{回転球殻 MHD ダイナモの定式化  1 基礎方程式系  3}

1 基礎方程式系

ここでは基礎方程式系の定式化について記述する. 流体はブシネスク近似を適応した磁気流体であるとし, 重力加速度は密度が一様な物体の自己重力場, 基本場の温度分布は内部一様熱源による熱伝導解であると仮定する. 速度場と磁場をトロイダル–ポロイダルポテンシャルへ分解し, その時間発展方程式と境界条件を定式化する (Zhang and Busse,1987, Chandrasekhar, 1961).

1.1 系の設定と支配方程式

角速度 Ω で回転する半径 r_{i}, r_{o} の球に挟まれた球殻内にある導電性流体について考える. 球殻内の流体はブシネスク流体かつ磁気流体であるとし, 座標系を球殻の中心を原点とする緯度経度座標系 (λ, φ, r) とする. これを図 1 に示す.

PSfrag replacements φ λ r r_{i} r_{o} Ω x y z

図 1: 考えている系の模式図支配方程式系は回転系におけるブシネスク流体の方程式系と磁気流体の方程式系とを組み合わせることで得られる. 以下にこれらを書き下す:

∇·u = 0,  (1)
∇·B = 0,  (2)

\frac{∂u}{∂t} + (u · ∇)u + 2Ωk × u = −\frac{1}{ρ_{0}} ∇p + \frac{ρ}{ρ_{0}} g + ν∇^{2}u + \frac{1}{μ_{0}ρ_{0}} (∇×B) × B,

(3)

\frac{∂T}{∂t} + (u · ∇)T = κ∇^{2}T + Q,  (4)
\frac{∂B}{∂t} = ∇×(u × B) + η∇^{2}B.  (5)

各々の記号が何を表しているかは次の通りである:

u  流速ベクトル.

成分毎に書く場合には,u = (u, v, w)

B 磁束密度ベクトル.

成分毎に書く場合には, B = (B_{λ}, B_{φ}, B_{r})

T  流体の実際の温度.

基本場の温度 T_{0}, そこからのずれ θ を用いて表すと, T = T_{0} + θ である.

ρ, ρ_{0} 流体の密度, 及びその基本場の密度.

ここではブシネスク流体であることを仮定しているので, ρ = ρ_{0}(1 − αθ) である.

p  流体の全圧力.

基本場の圧力 p_{0}, そこからのずれ π を用いて表すと p = p_{0} + π である.

k  球殻の回転軸方向の単位ベクトル ( 図 1 での z 軸方向).
g  重力加速度ベクトル.
Q  内部熱源 (温度上昇率).

また, 以下の物性定数は一定であるとする.

\begin{matrix}ν & 動粘性率. \\ κ & 熱拡散率. \\ α & 体積膨張率. \\ η & 磁気拡散率. \\ μ_{0} & 透磁率. \end{matrix} 重力加速度 g として密度がほぼ一様な物体の自己重力を考える事とし,

g = −γr  (6)

とする. ここで γ は定数, r は球殻の中心を原点とする動径方向の位置ベクトルである^{1}.

基本場は内部一様熱源分布 (Q = const) の下での磁場の無い静止状態であるとし, 基本場の温度分布として

T_{0}(r) = −\frac{1}{2}βr^{2} + C  (7)

を仮定する. だだし β ≡ \frac{1}{3} \frac{Q}{κ} とした^{2} .

この基本場を支配方程式系 (1)〜熱の式 (4) から分離すると,

∇·u = 0, ∇·B = 0  (8)

\frac{∂u}{∂t} + (u · ∇)u + 2Ωk × u = −\frac{1}{ρ_{0}} ∇π + αθγr + ν∇^{2}u + \frac{1}{μ_{0}ρ_{0}} (∇×B) × B,

(9)

\frac{∂θ}{∂t} + (u · ∇)θ − βu · r = κ∇^{2}θ,  (10)
\frac{∂B}{∂t} = ∇×(u × B) + η∇^{2}B  (11)

となる^{3} .

次にこれらの方程式系を以下のスケールを用いてスケーリングする:

長さスケール D = r_{o} − r_{i} (球殻の厚さ)
時間スケール D^{2}/ν  (粘性拡散時間)
速度スケール ν/D  (粘性速度)
温度スケール βD^{2}
磁場のスケール\sqrt{μ_{0}ρ_{0}}ν/D
圧力のスケール ρ_{0}ν^{2}/D^{2}

1重力加速度には基本場からの密度変化と回転の効果による非球対称な成分が存在するが, 簡単のためここでは無視する.

2熱伝導方程式を満たす温度分布には C′/r の不定性が残っているが, 後の数学的簡便さのため C′ = 0 の分布を選ぶこととする.

3基本場の温度分布として (7) を選んだことで, 基本場の移流が u · ∇T_{0} = −βu · r となる. すなわち, 本モデルでは基本場の温度分布がモデルに組み込まれている.

\overline{main.tex 平成 14 年 10 月 4 日 (佐々木洋平, 竹広真一, 林祥介)}

その結果, 方程式系は

∇·u = 0, ∇·B = 0,  (12)

\frac{∂u}{∂t} + (u · ∇)u + τk × u = −∇π + P_{r}^{−1}R_{a}θr + ∇^{2}u + (∇×B) × B,

(13)

\frac{∂θ}{∂t} + (u · ∇)θ − u · r = P_{r}^{−1}∇^{2}θ,  (14)
\frac{∂B}{∂t} − ∇×(u × B) = P_{m}^{−1}∇^{2}B  (15)

となる. ここで u, π, B, θ は, 各々無次元化された速度, 圧力のずれ, 磁場, 温度のずれを表すと定義しなおした.

スケーリングによって現れた無次元数は各々以下の通りである. R_{a} はレイリー数, τ はエクマン数 E の逆数, (もしくはテイラー数の平方根), P_{r} はプランドル数, P_{m} は磁気プランドル数であり, それぞれ次式で定義される^{4} :

R_{a} ≡ \frac{αγβD^{6}}{κν} ,  (16)
τ^{−1} = E ≡ \frac{ν}{2ΩD^{2}} ,  (17)
P_{r} ≡ \frac{ν}{κ},  (18)
P_{m} ≡ \frac{ν}{η} .  (19)

これに加えて系に現れる重要なパラメータとして内径外径比 ζ ≡ r_{i}/r_{o} が存在する.

4業界標準の無次元化と種々の論文での無次元化については, appendix 参照 (未整備).

\overline{main.tex 平成 14 年 10 月 4 日 (佐々木洋平, 竹広真一, 林祥介)}

1.2 ポテンシャルによる方程式系の表現

連続の式および磁場のソレノイダル条件 (12) より, 速度場と磁場をトロイダル, ポロイダルポテンシャルを用いて表すことができる:

u = ∇×(ψr) + ∇×∇×(ϕr),  (20)
B = ∇×(gr) + ∇×∇×(hr).  (21)

これらのポテンシャル ψ, ϕ, g, h を用いて方程式系を表す. 運動方程式 (13) に r·∇×, r · ∇×∇× を作用することで,

[(∇^{2} − \frac{∂}{∂t} ) L_{2} + τ \frac{∂}{∂λ} ] ψ − τQϕ = r · ∇×\{(u · ∇)u − (B · ∇)B)\} (22)

[(∇^{2} − \frac{∂}{∂t} ) L_{2} + τ \frac{∂}{∂λ} ] ∇^{2}ϕ + τQψ − P_{r}^{−1}R_{a}L_{2}θ

= −r · ∇×∇×{(u · ∇)u − (B · ∇)B)} (23)

が得られる. 式中の演算子は各々次式で定義した:

μ ≡ sin φ,  (24)
L_{2} ≡ −r^{2}∇^{2} + \frac{∂}{∂r} r^{2} \frac{∂}{∂r}
= −[\frac{1}{cos^{2} φ} \frac{∂^{2}}{∂λ^{2}} + \frac{1}{cos φ} \frac{∂}{∂φ} (cos φ)\frac{∂}{∂φ} ]
= −[\frac{1}{1 − μ^{2}} \frac{∂^{2}}{∂λ^{2}} + \frac{∂}{∂μ} \{(1 − μ^{2})\frac{∂}{∂μ} \}],  (25)
Q ≡ k · ∇ − \frac{1}{2}(L_{2}k · ∇ + k · ∇L_{2})
= −r sin φ∇^{2} + (L_{2} + r\frac{∂}{∂r} ) × (sin φ\frac{∂}{∂r} − \frac{cosφ}{r} \frac{∂}{∂φ} )
= −rμ∇^{2} + (L_{2} + r\frac{∂}{∂r} ) × (μ\frac{∂}{∂r} − \frac{1 − μ^{2}}{r} \frac{∂}{∂μ} ).  (26)

熱の式 (14) をポテンシャルで表すと

(P_{r}^{−1}∇^{2} − \frac{∂}{∂t} )θ + L_{2}ϕ = u · ∇θ.  (27)

誘導方程式のポテンシャルでの表現は (15) へ r · ∇×, r· を作用する事で得られる:

(P_{m}^{−1}∇^{2} − \frac{∂}{∂t} )L_{2}g = −r · ∇×∇×(u × B),  (28)
(P_{m}^{−1}∇^{2} − \frac{∂}{∂t} )L_{2}h = −r · ∇×(u × B).  (29)

\underline{回転球殻 MHD ダイナモの定式化  1 基礎方程式系  8}

r · ∇×, r · ∇×∇× の座標表現は

r · (∇×F) = \frac{1}{cos φ} \frac{∂F_{φ}}{∂λ} − \frac{1}{cos φ} \frac{∂}{∂φ} (F_{λ} cos φ),
r · (∇×∇×F) = \frac{1}{cos φ} \frac{∂(∇×F)_{φ}}{∂λ} − \frac{1}{cos φ} \frac{∂}{∂φ} [(∇×F)_{λ} cos φ]
= \frac{1}{cos φ} \frac{∂}{∂λ} (\frac{1}{r} \frac{∂}{∂r} (rF_{λ}) − \frac{1}{r cos φ} \frac{∂F_{r}}{∂λ} )

−\frac{1}{cos φ} \frac{∂}{∂φ} [cos φ (\frac{1}{r} \frac{∂F_{r}}{∂φ} − \frac{1}{r} \frac{∂}{∂r} (rF_{φ}))]

= \frac{1}{cos φ} \frac{∂}{∂λ} \frac{1}{r} \frac{∂}{∂r} (rF_{λ}) + \frac{1}{cos φ} \frac{∂}{∂φ} \frac{cos φ}{r} \frac{∂}{∂r} (rF_{φ})

−\frac{1}{r cos^{2} φ} \frac{∂^{2}F_{r}}{∂λ^{2}} − \frac{1}{r cos φ} \frac{∂}{∂φ} (cos φ\frac{∂F_{r}}{∂φ} )

= \frac{1}{r} \frac{∂}{∂r} [r (\frac{1}{cos φ} \frac{∂F_{λ}}{∂λ} + \frac{1}{cos φ} \frac{∂}{∂φ} (F_{φ} cos φ))] + \frac{L_{2}F_{r}}{r} .

1.2.1 速度場及び磁場とポテンシャルとの関係式トロイダル, ポロイダルポテンシャルから速度場を求める式は次の通りである:

u = \frac{∂ ψ}{∂φ} + \frac{1}{r cos φ} \frac{∂}{∂r} (r\frac{∂ ϕ}{∂λ}) ,  (30)
v = −\frac{1}{cos φ} \frac{∂ ψ}{∂λ} + \frac{1}{r} \frac{∂}{∂r} (r\frac{∂ ϕ}{∂φ}) ,  (31)
w = \frac{L_{2} ϕ}{r}.  (32)

同様に磁場の各成分について:

B_{λ} = \frac{∂ g}{∂φ} + \frac{1}{r cos φ} \frac{∂}{∂r} (r\frac{∂ h}{∂λ}) ,  (33)
B_{φ} = −\frac{1}{cos φ} \frac{∂ g}{∂λ} + \frac{1}{r} \frac{∂}{∂r} (r\frac{∂ h}{∂φ}) ,  (34)
B_{r} = \frac{L_{2} h}{r}.  (35)

逆に, 速度場の各成分からポテンシャルを求めるには, (20) に r·, r · ∇× を作用させればよい.

r · u = r · [∇×(ψr) + ∇×∇×(ϕr)] = r · ∇×∇×(ϕr) = L_{2}ϕ,

r · ∇×u = r · ∇×[∇×(ψr) + ∇×∇×(ϕr)] = r · ∇×∇×(ψr) = L_{2}ψ,

\overline{main.tex 平成 14 年 10 月 4 日 (佐々木洋平, 竹広真一, 林祥介)}

したがって,

ϕ = L_{2}^{−1}(rw), ψ = L_{2}^{−1}[r · (∇×u)].  (36)

ここで L_{2}^{−1} は L_{2} の逆演算子である. 同様に, (21) に r·, r · ∇× を作用させることにより磁場の各成分からポテンシャルを求める式が得られる.

h = L_{2}^{−1}(rB_{r}), g = L_{2}^{−1}[r · (∇×B)].  (37)

\underline{回転球殻 MHD ダイナモの定式化  1 基礎方程式系  10}

1.3 境界条件

1.3.1 運動学的境界条件運動学的境界条件は, 境界を横切る流れが存在しない, である. すなわち:

w = 0, at r = \frac{ξ}{1 − ξ} , \frac{1}{1 − ξ} .  (38)

これをポテンシャルで表現すれば,

ϕ = 0, at r = \frac{ξ}{1 − ξ} , \frac{1}{1 − ξ} .  (39)

1.3.2 力学的境界条件力学的境界条件としては, 応力無し (stress-free), もしくは滑べり無し (rigid) を考える. 応力無しの場合には

\frac{∂}{∂r} (\frac{u}{r}) = \frac{∂}{∂r} (\frac{v}{r} ) = 0, at r = \frac{ξ}{1 − ξ} , \frac{1}{1 − ξ} .  (40)

ポテンシャルで表現すれば

\frac{∂^{2}ϕ}{∂r^{2}} = \frac{∂}{∂r} (\frac{ψ}{r} ) = 0, at r = \frac{ξ}{1 − ξ} , \frac{1}{1 − ξ} .  (41)

滑べり無しの場合には,

u = v = 0, at r = \frac{ξ}{1 − ξ} , \frac{1}{1 − ξ} .  (42)

ポテンシャルで表現すれば

\frac{∂ϕ}{∂r} = ψ = 0, at r = \frac{ξ}{1 − ξ} , \frac{1}{1 − ξ} .  (43)

\overline{main.tex 平成 14 年 10 月 4 日 (佐々木洋平, 竹広真一, 林祥介)}

1.3.3 熱的境界条件

熱的境界条件としては, 境界において温度固定, もしくは熱フラックス固定を考える.

温度固定条件の場合には

θ = 0, at r = \frac{ξ}{1 − ξ} , \frac{1}{1 − ξ} .  (44)

熱フラックス固定の場合には^{5},

\frac{∂θ}{∂r} = 0, at r = \frac{ξ}{1 − ξ} , \frac{1}{1 − ξ} .  (45)

1.3.4 磁場の境界条件

磁場の境界条件として球殻の外部が不導体である場合と, 外部と内部の電気伝導度が等しい場合を考える.

球殻の外部が不導体の場合には磁場が境界において外部のポテンシャル場に接続される. すなわち, 外部の磁場を B_{e}, そのポテンシャルを W として

B = B_{e}, at r = \frac{ξ}{1 − ξ} , \frac{1}{1 − ξ} .  (46)

ただし B_{e} = ∇W, ∇^{2}W = 0 である.

球殻外部と内部の電気伝導度が等しい場合には, 磁場だけではなく電場の水平成分も境界において連続となる:

B = B_{e}, at r = \frac{ξ}{1 − ξ} , \frac{1}{1 − ξ} ,  (47)
E_{H} = E_{eH}, at r = \frac{ξ}{1 − ξ} , \frac{1}{1 − ξ}  (48)

ここで添字_{H} は水平成分を表している. 球殻外部の電磁場は別途与えられなければならない^{6}.

v_{i} = ω_{i} × r

5基本場の熱フラックスが境界において与えられている事になる.

6地球ダイナモで内核を解く際に良く用いられる方法は, 内核が角速度 ω_{i} = ω_{i}k で剛体回転しているとして, を誘導方程式 (15) へ代入し B_{e} を求める, である. 内核の角速度の時間変化は別途トルクの式を用いて求められる.

\overline{main.tex 平成 14 年 10 月 4 日 (佐々木洋平, 竹広真一, 林祥介)}

\underline{回転球殻 MHD ダイナモの定式化  2 水平離散化  12}

2 水平離散化

水平方向の離散化には擬スペクトル変換を用いる (Glatzmaier,1984). 線形項に関する水平微分については球面調和関数を基底とする直行関数展開によって評価し, 非線形項は格子点上で値を計算する.

2.1 スペクトル展開

水平方向の展開関数系として, 演算子 L_{2} の固有関数である球面調和関数 Y_{n}^{m}(λ, μ) を用いる. Y_{n}^{m}(λ, μ) はルジャンドル陪関数 P_{n}^{m}(μ) を用いて

Y_{n}^{m}(λ, μ) ≡ P_{n}^{m}(μ)e^{imλ}  (49)

で定義され, 演算子 L_{2} について

L_{2}Y_{n}^{m}(λ, μ) = n(n + 1)Y_{n}^{m}(λ, μ)  (50)

を満たす. ただし n ≥ |m| である.

ポテンシャル ψ の Y_{n}^{m}(λ, μ) による離散的な展開 ψ_{ij} は

ψ_{ij} ≡ ψ(λ_{i}, μ_{j}) = R_{e} ∑_{m=−N}^{N} ∑_{n=|m|}^{N} ψ̃_{n}^{m}Y_{n}^{m}(λ_{i}, μ_{j}),  (51)

その逆変換 ψ̃_{n}^{m} は

ψ̃_{n}^{m} = \frac{1}{4π} ∫_{−1}^{1} dμ ∫_{0}^{π} dλψ(λ, μ)Y_{n}^{m∗}(λ, μ)
= \frac{1}{I} ∑_{i=1}^{I} ∑_{j=1}^{J} ψ(λ_{i}, μ_{j})Y_{n}^{m∗}(λ_{i}, μ_{j})ω_{j}  (52)

である^{7}. 逆変換において積分を和に置き換えて評価する際は, λ 積分について Gauss の台形公式を, μ 積分については Gauss–Legendre の積分公式を用いた. μ_{j} は Gauss 緯度, ω_{j} は Gauss 荷重, λ_{i} は等間隔の格子である.

以下では簡単のため, 和の記号を ∑_{m=−N}^{N} ∑_{n=|m|}^{N} = ∑_{m,n}, ∑_{i=1}^{I} ∑_{j=1}^{J} = ∑_{i,j} と表す事とする.

7ここでは波数切断として三角形切断を用いている.

\overline{main.tex 平成 14 年 10 月 4 日 (佐々木洋平, 竹広真一, 林祥介)}

水平微分を含む項の格子点における値は, スペクトル展開により各々次式で求められる:

(\frac{∂ψ}{∂λ} )_{ij} = R_{e} ∑_{m,n} ∑_{i,j} imψ̃_{n}^{m}Y_{n}^{m}(λ_{i}, μ_{j}),  (53)
((1 − μ^{2})\frac{∂ψ}{∂μ} )_{ij} = R_{e} ∑_{m,n} ∑_{i,j} ψ̃_{n}^{m}(1 − μ^{2})\frac{∂}{∂μ} Y_{n}^{m}(λ_{i}, μ_{j}).  (54)

また, 非線形項を格子点上で評価する際に現れる水平微分を計算しておく:

( ˜\frac{1}{1 − μ^{2}} \frac{∂ψ}{∂λ} )_{n}^{m} = \frac{1}{4π} ∫_{−1}^{1} dμ ∫_{0}^{π} dλ(\frac{1}{1 − μ^{2}} \frac{∂ψ}{∂λ} )Y_{n}^{m∗}(λ, μ)
= \frac{1}{4π} ∫_{−1}^{1} dμ ∫_{0}^{π} dλ imψ(λ, μ)\frac{1}{1 − μ^{2}} Y_{n}^{m∗}(λ, μ)
= \frac{1}{I} ∑_{i,j} imψ_{ij}Y_{n}^{m∗}(λ_{i}, μ_{j})\frac{ω_{j}}{1 − μ^{2}}  (55)
( ˜\frac{∂ψ}{∂μ} )_{n}^{m} = \frac{1}{4π} ∫_{−1}^{1} dμ ∫_{0}^{π} dλ(\frac{∂ψ}{∂μ} )Y_{n}^{m∗}(λ, μ)
= −\frac{1}{4π} ∫_{−1}^{1} dμ ∫_{0}^{π} dλ ψ \frac{∂}{∂μ} Y_{n}^{m∗}(λ, μ)
= \frac{1}{I} ∑_{ij} ψ_{ij}(1 − μ_{j}^{2})\frac{∂}{∂μ} Y_{n}^{m∗}(λ_{i}, μ_{j})\frac{ω_{j}}{1 − μ^{2}} .  (56)

ここで Uψ が極 μ = ±1(φ = ±\frac{π}{2} ) において 0 となる事を用いている.

さらに L_{2} を含む項の評価には,

(L̃_{2}ψ)_{n}^{m} = n(n + 1)ψ̃_{n}^{m}  (57)

を用いる.

以上の計算結果は, ψ だけではなく ϕ, g, h, θ についてもそのまま成立する.

\underline{回転球殻 MHD ダイナモの定式化  2 水平離散化  14}

2.2 方程式のスペクトル表現

スペクトル展開されたポテンシャルで方程式系を書き表す. 2.2.1 運動方程式トロイダルポテンシャル:

\frac{∂}{∂t} ψ̃_{n}^{m} = G_{n}^{m}ψ̃_{n}^{m} − τ \frac{1}{n(n + 1)}(Q̃ϕ)_{n}^{m}
− \frac{1}{n(n + 1)}[r · ∇×(u · ∇u − B · ∇B)]_{n}^{m}. (n ≠ 0)  (58)

ポロイダルポテンシャル:

\frac{∂}{∂t} (D_{n}ϕ̃_{n}^{m}) = G_{n}^{m}(D_{n}ϕ̃_{n}^{m}) + τ \frac{1}{n(n + 1)}(Q̃ψ)_{n}^{m} − P_{r}^{−1}R_{a}θ̃_{n}^{m}
+ \frac{1}{n(n + 1)}[r · ∇×∇×(u · ∇u − B · ∇B)]_{n}^{m}. (n ≠ 0) (59)

ここで G_{n}^{m}, D_{n} は次式で定義される:

D_{n} ≡ \frac{∂}{∂r} + \frac{2}{r} \frac{∂}{∂r} − \frac{n(n + 1)}{r^{2}}, G_{n}^{m} ≡ \frac{im}{n(n + 1)} + D_{n}.  (60)

2.2.2 熱の式

\frac{∂}{∂t} θ̃_{n}^{m} = D_{n}θ̃_{n}^{m} + n(n + 1)θ̃_{n}^{m} − P_{r}^{−1}(u · ∇θ)_{n}^{m}  (61)

2.2.3 誘導方程式トロイダルポテンシャル

\frac{∂}{∂t} g̃_{n}^{m} = P_{m}^{−1}D_{n}g̃_{n}^{m} + \frac{1}{n(n + 1)}[r · ∇×(u × B)]_{n}^{m}. (n ≠ 0)  (62)

ポロイダルポテンシャル

\frac{∂}{∂t} h̃_{n}^{m} = P_{m}^{−1}D_{n}h̃_{n}^{m} − \frac{1}{n(n + 1)}[r · ∇×∇×(u × B)]_{n}^{m}. (n ≠ 0) (63)

\overline{main.tex 平成 14 年 10 月 4 日 (佐々木洋平, 竹広真一, 林祥介)}

2.2.4 速度場の表現

U_{ij}(r, t) = ∑_{m,n}\{ψ̃_{n}^{m}(1 − μ^{2})\frac{∂}{∂μ} Y_{n}^{m}(λ_{i}, μ_{j}) + im\frac{1}{r} \frac{∂}{∂r} (rϕ̃_{n}^{m})Y_{n}^{m}(λ_{i}, μ_{j})\}

(64)

V_{ij}(r, t) = ∑_{m,n}\{−imψ̃_{n}^{m}Y_{n}^{m}(λ_{i}, μ_{j}) + \frac{1}{r} \frac{∂}{∂r} (rϕ̃_{n}^{m})(1 − μ^{2})\frac{∂}{∂μ} Y_{n}^{m}(λ_{i}, μ_{j})\}

(65)

w_{ij}(r, t) = ∑_{m,n} \frac{n(n + 1)}{r}ϕ̃_{n}^{m}Y_{n}^{m}(λ_{i}, μ_{j})  (66)

2.2.5 磁場の表現

B_{λ_{ij}}(r, t) cos φ = ∑_{m,n}\{g̃_{n}^{m}(1 − μ^{2})\frac{∂}{∂μ} Y_{n}^{m}(λ_{i}, μ_{j}) + im\frac{1}{r} \frac{∂}{∂r} (rh̃_{n}^{m})Y_{n}^{m}(λ_{i}, μ_{j})\}

(67)

B_{μ_{ij}}(r, t) cos φ = ∑_{m,n}\{−img̃_{n}^{m}Y_{n}^{m}(λ_{i}, μ_{j}) + \frac{1}{r} \frac{∂}{∂r} (rh̃_{n}^{m})(1 − μ^{2})\frac{∂}{∂μ} Y_{n}^{m}(λ_{i}, μ_{j})\}

(68)

B_{r_{ij}}(r, t) = ∑_{m,n} \frac{n(n + 1)}{r}h̃_{n}^{m}Y_{n}^{m}(λ_{i}, μ_{j})  (69)

\underline{回転球殻 MHD ダイナモの定式化  3 鉛直離散化  16}

3 鉛直離散化

Glatzmaier(1984) に従い, 方程式を鉛直方向に離散化する. 鉛直微分はチェビシェフ関数展開を用いて評価する.

3.1 鉛直離散化

動径を差分化し, 上の層から下の層へと番号をつける.

格子点をチェビシェフポイント (チェビシェフ関数が極地をとる点) にとるため, 動径を次式のように座標変換する:

r = \frac{r_{o} + r_{i}}{2} + \frac{r_{o} − r_{i}}{2} ξ,  (70)
ξ = \frac{2r − (r_{o} + r_{i})}{r_{o} − r_{i}}.  (71)

この変換により, 動径は ξ を用いて区間 [−1, 1] で表現される (上部境界が 1, 下部境界が −1).

座標 ξ 上で K + 1 個の格子点を次のように定義する:

ξ_{k} = cos(\frac{πk}{K} ), k = 0, 1, ⋯ K.  (72)

各格子点でのチェビシェフ関数の値は

T_{l}(ξ_{k}) = cos(\frac{lkπ}{K} ).  (73)

格子点での物理量をチェビシェフ関数で展開する. 例えば ψ̃(r, t)_{n}^{m} について,

ψ_{nl}^{m}(t) = (\frac{2}{K} ) ∑_{k=0}^{K}′′ψ̃_{n}^{m}(ξ_{k}, t)T_{l}(ξ_{k}), l = 0, 1, ⋯ K.  (74)

ただし ∑′′ は最初と最後の項を 2 分の 1 倍して和をとることを表している. 同様にチェビシェフ展開係数から格子点での値を求めるには

ψ̃_{n}^{m}(ξ_{k}, t) = ∑_{l=0}^{K}′′ψ_{jn}^{m}(t)T_{j}(ξ_{k}).  (75)

以上の結果は, ϕ̃(r, t)_{n}^{m}, g̃(r, t)_{n}^{m}, h̃(r, t)_{n}^{m}θ̃(r, t)_{n}^{m} についてもそのまま成立する.

3.2 鉛直境界条件

内部領域の格子点において, 差分化した方程式を時間積分した後に, 境界条件を満たすように境界における格子点上の関数値を定める.

3.2.1 トロイダルポテンシャル

トロイダルポテンシャルに関しては r_{k} (k = 1, ⋯ K − 1 ) の点で時間積分を実行する.

応力無し条件の場合, 境界値は内部領域を計算した後

ψ̃_{n}^{m}(r_{0}, t) = \frac{r_{0}}{r_{1}} ψ̃_{n}^{m}(r_{1}, t), ψ̃_{n}^{m}(r_{K}, t) = \frac{r_{K}}{r_{K−1}} ψ̃_{n}^{m}(r_{K−1}, t)  (76)

で定める.

滑べり無し条件の場合は

ψ̃_{n}^{m}(r_{0}, t) = 0, ψ̃_{n}^{m}(r_{K}, t) = 0.  (77)

となる.

3.2.2 ポロイダルポテンシャル

ポロイダルポテンシャルは r_{k}, (k = 2, ⋯ K − 2) の点で (D_{n}ϕ̃_{n}^{m}) について時間積分を実行した後, D_{n} の逆演算 (積分) をおこなって ϕ̃_{n}^{m} を定める. そのためには境界値問題

(\frac{d^{2}}{dr^{2}} + \frac{2}{r} \frac{d}{dr} − \frac{n(n + 1)}{r^{2}})ϕ̃_{n}^{m} = (D_{n}ϕ̃_{n}^{m}),  (78)
ϕ̃_{n}^{m} = \frac{∂^{2}ϕ̃_{n}^{m}}{∂r^{2}} = 0 at r = r_{i}, r_{o}.  (79)

を解かねばならない (上式は応力無し条件の場合). 滑べり無し条件の場合には、 (79) を

ϕ̃_{n}^{m} = \frac{∂ϕ̃_{n}^{m}}{∂r} = 0 at r = r_{i}, r_{o}.  (80)

へ置き換えればよい.

\underline{回転球殻 MHD ダイナモの定式化  3 鉛直離散化  18}

(78), (79) をチェビシェフ関数展開で表現すると

∑_{l=0}^{K}′′(\frac{d^{2}T_{j}(r_{k})}{dr^{2}} + \frac{2}{r_{k}} \frac{dT_{j}(r_{k})}{dr} − \frac{n(n + 1)}{r_{k}^{2}}T_{j}(r_{k}))ϕ_{jn}^{m} = D_{n}ϕ̃_{n}^{m}(r_{k}),

k = 2, ⋯ K − 2,  (81)
∑_{l=0}^{K}′′T_{j}(r_{0})ϕ_{jn}^{m} = 0, ∑_{l=0}^{K}′′T_{j}(r_{K})ϕ_{jn}^{m} = 0,  (82)
∑_{l=0}^{K}′′ \frac{d^{2}T_{j}(r_{0})}{dr^{2}}ϕ_{jn}^{m} = 0, ∑_{l=0}^{K}′′ \frac{d^{2}T_{j}(r_{K})}{dr^{2}}ϕ_{jn}^{m} = 0.  (83)

と, ϕ_{jn}^{m} に関する連立 1 次方程式となる. これを LU 分解で解き, ϕ_{jn}^{m} を定める. 滑べり無し条件の場合は (83) が

∑_{l=0}^{K}′′ \frac{dT_{j}(r_{0})}{dr}ϕ_{jn}^{m} = 0, ∑_{l=0}^{K}′′ \frac{dT_{j}(r_{K})}{dr}ϕ_{jn}^{m} = 0,  (84)

となる. 3.2.3 温度温度に関しては r_{k}, k = 1, ⋯ K − 1 の点で時間積分を実行する. 境界値は内部領域を計算した後 (温度固定の場合)

Θ̃_{n}^{m}(r_{0}, t) = 0, Θ̃_{n}^{m}(r_{K}, t) = 0,  (85)

で定める. 熱フラックス固定の場合は (85) を

Θ̃_{n}^{m}(r_{0}, t) = Θ̃_{n}^{m}(r_{1}, t), Θ̃_{n}^{m}(r_{K}, t) = Θ̃_{n}^{m}(r_{K−1}, t),  (86)

に置き換える. 3.2.4 磁場先ず境界外部が不導体の場合を考える. 接続すべきポテンシャル場 W を球面調和函数展開すると,

W = ∑_{n,m} W_{n}^{m}r^{−n−1}Y_{n}^{m} (r > r_{o}),  (87)
= ∑_{n,m} W_{n}^{m}r^{n}Y_{n}^{m} (r < r_{i}).  (88)

\overline{main.tex 平成 14 年 10 月 4 日 (佐々木洋平, 竹広真一, 林祥介)}

従って各成分は

B_{λ} cos θ = ∑_{n,m} imW_{n}^{m}r^{−n−2}Y_{n}^{m} (r > r_{o}),  (89)
= ∑_{n,m} imW_{n}^{m}r^{n−1}Y_{n}^{m} (r < r_{i}),  (90)
B_{μ} cos θ = ∑_{n,m} W_{n}^{m}r^{−n−2}(1 − μ^{2})\frac{∂Y_{n}^{m}}{∂μ} (r > r_{o}),  (91)
= ∑_{n,m} W_{n}^{m}r^{n−1}(1 − μ^{2})\frac{∂Y_{n}^{m}}{∂μ} (r < r_{i}),  (92)
B_{r} = − ∑_{n,m}(n + 1)W_{n}^{m}r^{−n−2}Y_{n}^{m} (r > r_{o}),  (93)
= ∑_{n,m} nW_{n}^{m}r^{n−1}Y_{n}^{m} (r < r_{i}).  (94)

これに対して, 磁場の各成分が連続でなければならない. 一方, 磁場の鉛直成分 (69) と比較することにより,

\frac{n(n + 1)h_{n}^{m}}{r} = −(n + 1)W_{n}^{m}r^{−n−2} (r = r_{o}),  (95)
\frac{n(n + 1)h_{n}^{m}}{r} = nW_{n}^{m}r^{n−1} (r = r_{i}).  (96)

さらに水平成分を (67), (68) と比較して

g_{n}^{m} = 0, \frac{1}{r} \frac{∂(rh_{n}^{m})}{∂r} = W_{n}^{m}r^{−n−2}, (r = r_{o}),  (97)
g_{n}^{m} = 0, \frac{1}{r} \frac{∂(rh_{n}^{m})}{∂r} = W_{n}^{m}r^{n−1}, (r = r_{i}),  (98)

W_{n}^{m} を代入して, h_{n}^{m} の境界条件が得られる:

\frac{∂h_{n}^{m}}{∂r} + \frac{(n + 1)h_{n}^{m}}{r} = 0, (r = r_{o}),  (99)
\frac{∂h_{n}^{m}}{∂r} − \frac{nh_{n}^{m}}{r} = 0, (r = r_{i}),  (100)

従ってまとめると, トロイダル磁場に関しては r_{k}, k = 1, ⋯ K − 1 の点で時間積分を実行する. 境界値は内部領域を計算した後,

g̃_{n}^{m}(r_{0}, t) = 0, g̃_{n}^{m}(r_{K}, t) = 0,  (101)

で定める. ポロイダル磁場も r_{k}, k = 1, ⋯ K − 1 の点で時間積分を実行する. その後境界条件により境界の値を定める:

\frac{∂h_{n}^{m}}{∂r} + \frac{(n + 1)h_{n}^{m}}{r} = 0, at r = r_{0},  (102)
\frac{∂h_{n}^{m}}{∂r} − \frac{nh_{n}^{m}}{r} = 0, at r = r_{K}.  (103)

また, この時の外部磁場は次式で表される:

W_{n}^{m} = −n(r^{n+1}h_{n}^{m})_{r=r_{0}} (r > r_{o}), W_{n}^{m} = [\frac{(n + 1)h_{n}^{m}}{r^{n}} ]_{r=r_{K}}, (r < r_{i}). (104)

次に球殻内の流体と球殻の芯となる領域 (以下, 内核) との電気伝導度が等しい場合を, 特に内核が系に対して角速度 ω^{i} = ω^{i}k で剛体回転している場合について考える (以下, 添字^{i} は内核を表す).

内核の磁場 B^{i} 及び角速度 ω_{i} は, 剛体回転する導体の誘導方程式と自転軸周りのトルクの釣合から,

\frac{∂B^{i}}{∂t} = ω^{i} \frac{∂B^{i}}{∂λ} + P_{m}^{−1}∇^{2}B^{i}  (105)
I \frac{∂ω^{i}}{∂t} = Γ_{M} + Γ_{ν}  (106)

によって求められる. ここで I は内核の自転軸周りの慣性モーメントであり, Γ_{M} はローレンツ力によるトルクを, Γ_{ν} は粘性応力によるトルクを表す.

• 回転対称性を用いる:

B^{i}(λ, φ, −r) = B^{i}(λ ± π, π − φ, r)

• h_{n}^{mi}(−r) = (−1)^{l+1}h_{n}^{mi}(r)

• チェボシェフで計算するのは, r = [0, r] まで.

• 境界条件は rigid →ポテンシャル, ポテンシャルの動径微分が連続

• トルクは…むーん.

\underline{回転球殻 MHD ダイナモの定式化参考文献 21}

参考文献

Chandrasekhar,S.,1961: Hydrodynamic and Hydromagnetic Stability. Oxford University Press.

Glatzmaier,G.A.,1984: Numerical simulations of stellar convective dynamos. I.the model and method. J. Comput. Phys., 55, 461–484.

Hollerbach,R.,2000: A spectral solution of the magneto-convection equations in spherical geometry. J. Numer. Mech. Fluids, 32, 773-797.

Zhang,K. and Busse,F.H.,1987: On the onset of convection in rotating spherical shells. Geophys. Astrophys. Fluid Dyn., 39, 119–147.