%%% 表題: DM2ゼミ 発表資料 %%% %%% 履歴: 2004/05/06 杉山耕一朗 %%% 履歴: 2004/05/09 杉山耕一朗 %%% %include "default.mgp" %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %nodefault %fore "red", size 8, back "darkblue", font "thick", vgap 0, ccolor "gray" %bgrad 0 0 128 0 1 "white" "white" "blue" "white" "white" "white" "white" "white" %center, fore "yellow", font "thick" 木星型惑星大気の熱力学計算 %size 6 最適化法によるギブス自由エネルギーの最小化 %size 5, fore "black", font "thick" 杉山耕一朗(北大理) sugiyama@gfd-dennou.org 2004 年 6 月 16 日 於 DM2 ゼミ %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 目次 木星大気の大気構造 平衡雲凝縮モデルの紹介 Weidenschilling and Lewis (1973), Atreya and Romani (1985) 地球流体電脳倶楽部版平衡熱力学計算コード oboro 想定する系 計算手順 化学平衡計算手法 歴史 最適化法の紹介 RAND 法の紹介 RAND 法を化学平衡計算に適用するにあたって 計算結果の紹介 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" 木星の大気構造 光学観測 表層のアンモニア雲層のみ観測可能 深部の組成, 温度構造はわからない 電波観測 得られた輝度温度と一致するような大気構造を推定 直接的に深部の組成, 温度構造が得られるわけでない ガリレオプローブの観測 雲のない乾燥領域に落下 雲に覆われた湿潤領域の情報は得られてない %fore "red", center, size 6 大気構造を数値モデルによって考察 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" 平衡雲凝縮モデル (ECCM) %size 2 断熱線に沿った鉛直 1 次元の %cont, fore "red" 温度構造, 組成変化 %cont, fore "black" の計算 Weidenschilling and Lewis (1973) Atreya and Romani (1985) 雲分布の %cont, fore "red" 「標準推定値」 %fore "black" %center %image "images/WL1973f3.ps" 0 100 100 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" ECCM の定式化 断熱線の計算 dX, v を dp, dT で書き換え, dT/dz の式を解く dp は静水圧平衡より dz へ変換 %% %size 2 %filter "platex2eps.sh eccm" \color[named]{Red} {\boldmath \begin{eqnarray} dS = \bar{c}_{p} dT - v dp + \sum L_{k} dX_{k} + \sum L_{r} dX_{r} = 0 \nonumber \end{eqnarray} } %endfilter %mark %center, image "eccm.eps" 0 250 250 1 %left 不便な点 大気中で生じ得る化学反応を知らねば解けない 水溶液の扱いが面倒 変数として水溶液濃度が追加 濃度を仮定して飽和蒸気圧を計算. 飽和蒸気圧から予想される凝縮量と仮定した濃度が 等しくなるまで反復計算 %filter "platex2eps.sh eccm_setumei" \color[named]{Red} {\boldmath ~\\ $T$: 温度\\ $p$: 圧力\\ $\bar{c_{p}}$: 平均比熱\\ $v$: 体積\\ $L_{k}$: 潜熱\\ $dX_{k}$: 凝結によるモル変化\\ $L_{r}$: 反応熱\\ $dX_{r}$: 反応によるモル変化\\ } %endfilter %again %right, image "eccm_setumei.eps" 0 150 150 1 %cont %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" 研究目的 %center, fore "red", size 6 木星型惑星大気を想定した 平衡熱力学計算コードの開発 %left, fore "black", size 3 計算コード開発にあたって 多様な惑星大気に対して応用可能 様々な温度・圧力条件で計算可能 複数の化学種の相変化, 化学反応, 混合を簡単に扱えるように ECCM に代わるものを開発 ギブス自由エネルギー最小化法を採用 静的安定度 N^2 を見積もる 惑星大気のキーパラメタの 1 つ 木星型惑星大気の静的安定度も湿潤対流によって決まると仮定 N: 浮力振動数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" 我々の計算 擬断熱変化を仮定 平衡雲凝縮モデル(ECCM)と同様の仮定 計算方法は ECCM とは異なる %fore "darkgreen", size 5 --> %cont, fore "red" ギブス自由エネルギー最小化法 %fore "black" 計算手順 温度圧力を固定し, 元素数保存の条件のもとギブス自由エネルギーが最小化となる平衡組成を求める. 当該温度圧力における平衡組成を用いてエントロピーを計算する. 断熱線 dS = 0 をエントロピーの保存される温度圧力を探索することで求める. 断熱線 dS = 0 に沿った種々の大気科学的な物理量の計算 静的安定度 N^2 の見積もり %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" 平衡組成の計算 与えられた温度, 圧力, 元素組成に対して, %cont, fore "red" ギブス自由エネルギーを最小化 %fore "black" 大気中で生じ得る %cont, fore "red" 化学反応 %cont, fore "black" の詳細を知らなくて済み, 元素組成を変化させても数値計算コードを書きなおす必要がほとんどない. ギブス自由エネルギーの定義 %filter "platex2eps.sh gibbs" \color[named]{Red} {\boldmath \begin{eqnarray} G(T, p, n^{\phi}_{i}) &=& \sum \mu_{i}^{\phi}(T, p, n^{\phi}_{i}) n_{i}^{\phi} \nonumber \\ &=& \sum \left\{ {\mu_{i}^{\circ}}^{\phi}(T) + RT \ln \frac{n_{i}^{\phi}}{\sum n_{i}^{\phi}} + \alpha RT \ln{\frac{p}{p_0}} \right\} n_{i}^{\phi} \nonumber \end{eqnarray} } %endfilter %center %image "gibbs.eps" 0 250 250 1 %left %mark %fore "black" パラメタ 温度 圧力 基準状態での化学ポテンシャル %again %filter "platex2eps.sh gibbs_hensuu1" \color[named]{Red} {\bf \boldmath ~\\ $T$: temperature \\ $p$: pressure \\ $\mu$: chemical potential \\ $n$: composition \\ $R$: gas constant } %endfilter %filter "platex2eps.sh gibbs_hensuu2" \color[named]{Red} {\bf \boldmath ~\\ $^{\circ}$: standard state \\ $\phi$: phase \\ \\ if solution or mixing gas, $\alpha = 1$ \\ otherwise $\alpha = 0$ } %endfilter %right %image "gibbs_hensuu1.eps" 0 150 150 1 %cont %image "gibbs_hensuu2.eps" 0 150 150 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" エントロピーの計算 平衡組成に対するエントロピー %filter "platex2eps.sh entropy" \color[named]{Red} {\bf \boldmath \begin{eqnarray} S = - \left( \DP{G}{T} \right)_{p, n_{i}} = - \sum_{i} \left\{ \DP{{\mu_i^{\circ}}^{\phi}(T)}{T} + R \ln \left( \frac{n_i^{\phi}}{\sum n_i^{\phi}}\right) + \alpha R \ln \frac{p}{p_{0}} \right\} n_{i}^{\phi} \nonumber \end{eqnarray} } %endfilter %center %image "entropy.eps" 0 250 250 1 %left %fore "black" パラメタ 温度 圧力 基準状態での化学ポテンシャル 平衡組成 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" 擬断熱変化 擬断熱線 dS = 0 エントロピーが保存するような温度と圧力の関係を求める %center %image "images/pseudo3-1.eps" 0 180 180 1 %cont %image "images/pseudo3-2.eps" 0 180 180 1 %cont %image "images/pseudo3-3.eps" 0 180 180 1 %cont %image "images/pseudo3-4.eps" 0 180 180 1 %left %size 3 Step 1: 初期温度での平衡組成を計算し, エントロピー $S_0$ を求める. Step 2: 圧力を $dp$ だけ変化させる. 温度を変化させた時のエントロピーを 順次計算し, 前のステップでのエントロピー $S_0$ と一致する温度を求める. Step 3: ($T_{1}, p_{1}$) において凝縮が生じた場合, 次のステップにおいて 保存させるエントロピーは ($T_{1}, p_{1}$) での気相のエントロピー Step 4: Step 1 から Step 3 において得られた温度・圧力を結んで 断熱線を引く. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" 大気科学的な物理量の計算 静的安定度 乾燥断熱減率と乾燥断熱減率との差から評価 地球では, 空気塊は下降する際に放射冷却を受け, 大気全体の温度構造は湿潤断熱的になる そのような環境場で乾燥空気塊の安定性を調べる %filter "platex2eps.sh N2" \color[named]{Red} {\bf \boldmath \begin{eqnarray} N^{2} = \frac{g}{T_{v}}\left(\Gamma_{\rm m} - \Gamma_{\rm d}\right) \nonumber \end{eqnarray} } %endfilter %filter "platex2eps.sh N2-hosoku" \color[named]{Red} {\bf \boldmath ~\\ $T_{v}$: virtual temperature \\ $g$: gravity of planet \\ $\Gamma_{\rm m}$: adiabatic lapse rate of moist air \\ $\Gamma_{\rm d}$: adiabatic lapse rate of dry air } %endfilter %center %image "N2.eps" 0 250 250 1 %cont %cont %image "N2-hosoku.eps" 0 200 200 1 %image "images/cloud.eps" 0 150 150 1 %size 3 地球での温度構造の決まり方 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %nodefault %fore "red", size 8, back "darkblue", font "thick", vgap 0, ccolor "gray" %bgrad 0 0 128 0 1 "white" "white" "blue" "white" "white" "white" "white" "white" %center, fore "yellow", font "thick" 化学平衡計算手法の紹介 %size 5 最適化法によるギブス自由エネルギー最小化: RAND 法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" はじめに ギブス自由エネルギー最小化法 変数として, 温度, 圧力, 組成を選択 平衡状態をギブス自由エネルギーの最小化された状態とする 温度・圧力を与え, 元素数保存の条件の下でギブス自由エネルギーを最小化 %% %filter "platex2eps.sh ming" \color[named]{Red} {\boldmath \begin{eqnarray} && G(T, p, \Dvect{n}) = \sum_{i=1}^{\rm N} \sum_{\phi=1}^{\rm p} n_{i}^{\phi} \mu_{i}^{\phi}(T, p, \Dvect{n}) \nonumber \\ &&\sum_{\phi}\sum_{i} a_{ie}^{\phi} n_{i}^{\phi} = B_{e} \nonumber \end{eqnarray} } %endfilter %% %filter "platex2eps.sh ming-2" \color[named]{Red} {\boldmath $a_{ie}$: 化学種 $i$ に含まれる元素 $e$ の個数, $B_{e}$: 元素 $e$ の総量 } %endfilter %center %image "ming.eps" 0 250 250 1 %image "ming-2.eps" 0 200 200 1 %left oboro では 最小化の数値計算手法として RAND 法を利用 「条件つき最適化法」の「Newton 法」に分類 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 化学平衡計算 平衡問題の困難 方程式が非線形 熱力学関数の極値が 1 つとは限らない 局所的な極値も存在 進化を支えたのは 計算機の能力の発展 応用数学的な最適化問題の数値的取り扱い方法の確立 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 化学工学論文での計算手法の変遷 60 年代まで 「最適化法(optimize method)」 ギブス自由エネルギーの最小化される平衡組成を探索 様々な最適化法が化学平衡問題に応用されていた 「非線形方程式の解(solution of non-linear equation)」 相の間で化学ポテンシャル一定の条件を解く 70 -- 80 年代 「最適化法」がよく利用されるように 相の安定性解析手法の研究 80 年代後半 -- 現在 「最適化法」の「準 Newton 法」の「BFGS 公式」が標準 local minimum を避けるために確率的要素も加味 モンテカルロシミュレーション minimum が得られたら, 確率的に組成を振る %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 最適化法 各種の制約下で目的関数の極値を求める手法 制約なし最適化法 制約つき最適化法 最適化アルゴリズムの一般形 初期近似解 x_0 から計算開始 反復公式によって最適解に収束する点列 x_k を求める 「探索方向」, 「ステップ幅」の選び方でアルゴリズムが決まる %filter "platex2eps.sh optimize" \color[named]{Red} {\boldmath \begin{eqnarray} && \Dvect{x}_{k+1} = \Dvect{x}_{k} + \alpha_{k} \Dvect{d}_{k} \nonumber \end{eqnarray} } %endfilter %filter "platex2eps.sh optimize_setumei" \color[named]{Red} {\boldmath $\alpha_{k}$: ステップ幅, $\Dvect{d}_{k}$: 探索方向 } %endfilter %center %image "optimize.eps" 0 350 350 1 %image "optimize_setumei.eps" 0 200 200 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" Newton 法 アルゴリズム ステップ幅は 1 とする 探索方向は以下の連立一次方程式より求める %center %filter "platex2eps.sh newton" \color[named]{Red} {\boldmath \begin{eqnarray} Q &=& f(\Dvect{x}_{k}) + \Dgrad f(\Dvect{x}_{k}) (\Dvect{x} - \Dvect{x}_{k}) + \Dinv{2} \Dlapla f(\Dvect{x}_{k}) (\Dvect{x} - \Dvect{x}_{k})^{2} \nonumber \\ \DD{Q}{\Dvect{x}} &=& \DD{}{\Dvect{x}} f(\Dvect{x}_{k}) + \Dgrad f(\Dvect{x}_{k}) \cdot \DD{}{\Dvect{x}} (\Dvect{x} - \Dvect{x}_{k}) + \Dinv{2} \Dlapla f(\Dvect{x}_{k}) \cdot \DD{}{\Dvect{x}} (\Dvect{x} - \Dvect{x}_{k})^{2} \nonumber \\ &=& \Dgrad f(\Dvect{x}_{k}) + \Dgrad^{2} f(\Dvect{x}_{k}) \cdot (\Dvect{x} - \Dvect{x}_{k}) \nonumber \end{eqnarray} } %endfilter %% %filter "platex2eps.sh newton-H" \color[named]{Red} {\boldmath \begin{eqnarray} \Dgrad f (\Dvect{x}_{k}) + \Dlapla f(\Dvect{x}_{k}) \Dvect{d}_{k} = 0 \nonumber \end{eqnarray} } %endfilter %image "newton-H.eps" 0 300 300 1 %left 導出 目的関数を近似解 x_k のまわりでテーラー展開し, 2 次の微小量まで考慮 %image "newton.eps" 0 250 250 1 %filter "platex2eps.sh lapla" {\boldmath \begin{eqnarray} \Dlapla f(\Dvect{x}_{k}) \nonumber \end{eqnarray} } %endfilter %filter "platex2eps.sh grad" {\boldmath \begin{eqnarray} \Dgrad f(\Dvect{x}_{k}) \nonumber \end{eqnarray} } %endfilter %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" Newton 法 探索ベクトルのイメージ %center %image "images/optimize.eps" 0 250 250 1 %size 3 適当な組成のまわりで 2 次関数に近似しその極値を反復的に 求めることで, ポテンシャルの最も低い状態へ収束させる %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 準ニュートン法 Newton 法の探索ベクトル中の %cont %image "lapla.eps" 0 250 250 1 %cont を直接計算しない方法の総称 x_k と %cont %image "grad.eps" 0 250 250 1 %cont の情報を用いて %cont %image "lapla.eps" 0 250 250 1 %cont 近似行列を作成 BFGS 公式 現在, 経験的にも理論的にも最も優れているとされている準 Newton 法 %center %filter "platex2eps.sh bfgs" \color[named]{Red} {\boldmath \begin{eqnarray} H_{k+1} = H_{k} - \frac {H_{k} \Dvect{s}_{k} \Dvect{s}_{k}^{T} H_{k}} {\Dvect{s}_{k}^{T} H_{k} \Dvect{s}_{k} } + \frac {\Dvect{y}_{k} \Dvect{y}_{k}^{T}} {\Dvect{s}_{k}^{T} \Dvect{y}_{k} } \nonumber \end{eqnarray} } %endfilter %filter "platex2eps.sh hessian" \color[named]{Red} {\boldmath \begin{eqnarray} \Dvect{d}_{k} \equiv - H_{k}^{-1} \Dgrad f(\Dvect{x}) \nonumber \end{eqnarray} } %endfilter %filter "platex2eps.sh s_y" \color[named]{Red} {\boldmath \begin{eqnarray} \Dvect{s}_{k} \equiv \Dvect{x}_{k+1} - \Dvect{x}_{k}, \;\;\; \Dvect{y}_{k} \equiv \Dgrad f(\Dvect{x}_{k+1}) - \Dgrad f(\Dvect{x}_{k}) \nonumber \end{eqnarray} } %endfilter %image "hessian.eps" 0 250 250 1 %image "bfgs.eps" 0 250 250 1 %image "s_y.eps" 0 200 200 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 制約つき Newton 法 %size 3 アルゴリズム ステップ幅は 1 探索方向は制約なし Newton 法と同様に決める %% %filter "platex2eps.sh lambda" \color[named]{Red} {\boldmath \begin{eqnarray} && \Dgrad_{x} L(\Dvect{x}, \Dvect{\lambda}) = \Dgrad_{x} f(\Dvect{x}) - (J(\Dvect{x}))^{T} \Dvect{\lambda} = 0 \nonumber \\ && \Dgrad_{\lambda} L(\Dvect{x}, \Dvect{\lambda}) = \Dvect{g}(\Dvect{x}) = 0 \nonumber \end{eqnarray} } %endfilter %% %% %filter "platex2eps.sh newton_con" \color[named]{Red} {\boldmath \[ \left[ \begin{array}{cc} H(\Dvect{x}_{k}, \Dvect{\lambda}_{k}) & -J^{T}(\Dvect{x}_{k}) \\ J(\Dvect{x}_{k}) & 0 \end{array} \right] \left[ \begin{array}{c} \Dvect{d}_{k} \\ \Dvect{\lambda}_{k+1} \end{array} \right] = - \left[ \begin{array}{c} \Dgrad_{x} f(\Dvect{x}_{k}) \\ \Dvect{g}(\Dvect{x}_{k}) \end{array} \right] \] } %endfilter %filter "platex2eps.sh jacobi" \color[named]{Red} {\boldmath \begin{eqnarray} J(\Dvect{x}) \equiv \Dgrad \Dvect{g}(\Dvect{x}) \nonumber \end{eqnarray} } %endfilter %% %center %image "newton_con.eps" 0 250 250 1 %image "jacobi.eps" 0 250 250 1 %left 導出 条件つき極値を求める ラグランジュの未定乗数法を利用 %center %image "lambda.eps" 0 250 250 1 %left 特徴 収束域が狭い 初期値 x_0 の選び方によっては発散したり, 局所的な極値に陥りやすい %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 我々の計算: RAND 法 RAND 法 ギブス自由エネルギーを最小化するための手法として開発 White et al (1958) 条件: 元素数保存 「条件つき最適化法」の「Newton 法」に分類 Newton 法との差異 目的関数がギブス自由エネルギー自体 解くべき連立方程式の変数が, 探索ベクトルではなく, 最適解自体 等式制約が組成の一次関数なので, 等式制約をそのまま使える最適解の連立方程式として整理する方が都合が良い. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" RAND 法の定式化 最適解の決め方 条件式: Ax = B の形式 %% %filter "platex2eps.sh rand-1" \color[named]{Red} {\boldmath \begin{eqnarray} &&Q(\Dvect{n}_{k+1}) = G(\Dvect{n}_{k}) + \Dgrad{G(\Dvect{n}_{k}) } \cdot (\Dvect{n}_{k+1} - \Dvect{n}_{k}) + \Dinv{2} \Dlapla{G(\Dvect{n}_{k}) } \cdot (\Dvect{n}_{i+1} - \Dvect{n}_{k})^{2}. \nonumber \\ &&L = Q - \Dvect{\lambda} \left( \Dvect{A} \Dvect{n}_{k+1} - \Dvect{B}\right) = \Dgrad{G(\Dvect{n}_{k}) } + \Dlapla{G(\Dvect{n}_{k}) } \cdot (\Dvect{n}_{i+1} - \Dvect{n}_{k}) - \Dvect{\lambda} \left( \Dvect{A} \Dvect{n}_{k+1} - \Dvect{B}\right) \nonumber \\ &&\DP{Q(\Dvect{n}_{k+1})}{\Dvect{n}_{k+1}} = \left\{ \Dgrad{G(\Dvect{n}_{k})} - \Dlapla{G(\Dvect{n}_{k}) } \Dvect{n}_{k} \right\} + \Dlapla{G(\Dvect{n}_{k}) } \Dvect{n}_{i+1} - \Dvect{\lambda} \Dvect{A} = 0 \nonumber \end{eqnarray} } %endfilter %% %filter "platex2eps.sh rand-2" \color[named]{Red} {\boldmath \[ \left[ \begin{array}{cc} \Dlapla{G(\Dvect{x}_{k})} & \Dvect{A} \\ \Dvect{A} & 0 \end{array} \right] \left[ \begin{array}{c} \Dvect{n}_{k+1}\\ \Dvect{\lambda} \end{array} \right] = \left[ \begin{array}{c} \Dgrad{G(\Dvect{n}_{k})} - \Dlapla{G(\Dvect{n}_{k}) } \Dvect{n}_{k}\\ \Dvect{B} \end{array} \right] \] } %endfilter %center %image "rand-2.eps" 0 250 250 1 %left 導出 ギブス自由エネルギーを適当な組成のまわりで展開 %center %image "rand-1.eps" 0 200 200 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" RAND 法の定式化 方程式数を減らす 元素保存式をうまく利用 方程式数: 相の数 x 化学種の数 + 元素数 --> 元素数 + 相の数 %% %filter "platex2eps.sh rand-3" \color[named]{Red} {\boldmath % \[ \left[ \begin{array}{ccccc} \sum_{i} a_{i1} a_{i1} m_{i} & \sum_{i} a_{i2} a_{i1} m_{i} & \cdots & \sum_{i} a_{ie} a_{i1} m_{i} & B_{1} \\ \sum_{i} a_{i1} a_{i2} m_{i} & \sum_{i} a_{i2} a_{i2} m_{i} & \cdots & \sum_{i} a_{ie} a_{i2} m_{i} & B_{2} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ \sum_{i} a_{i1} a_{if} m_{i} & \sum_{i} a_{i2} a_{if} m_{i} & \cdots & \sum_{i} a_{ie} a_{if} m_{i} & B_{f} \\ B_{1} & B_{2} & \cdots & B_{e} & 0 \end{array} \right] \left[ \begin{array}{c} \lambda_{1} \\ \lambda_{2} \\ \vdots \\ \lambda_{e} \\ \left(\frac{\sum n_{i}}{\sum m_{i}} - 1 \right) \end{array} \right] = \left[ \begin{array}{c} \sum a_{i1} m_{i} \frac{\mu_{i}}{RT} \\ \sum a_{i2} m_{i} \frac{\mu_{i}}{RT} \\ \vdots \\ \sum a_{if} m_{i} \frac{\mu_{i}}{RT} \\ \sum m_{i} \frac{\mu_{i}}{RT} \\ \end{array} \right] \] % } %endfilter %% %filter "platex2eps.sh rand-4" \color[named]{Red} {\boldmath \begin{eqnarray} n_{i} = m_{i} \left\{ - \frac{\mu_{i}}{RT} + \frac{\sum n_{i}}{\sum m_{i}} + \sum_{e} \lambda_{e} a_{ie} \right\}. \nonumber \end{eqnarray} } %endfilter %center %image "rand-4.eps" 0 250 250 1 %image "rand-3.eps" 0 200 200 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 計算がうまく動かない! 局所的な極値へ収束 %center, image "images/H2O_on.gif" 0 100 100 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 計算がうまく動かない! エントロピーが収束しない %image "./images/gibbs.eps" 0 200 200 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 計算がうまく動かない! すぐに係数行列が非正則になる %size 6 ---> 計算自体が落ちるので絵無し %pause %center, fore "red" 最適化法を適用する系に則したチューニングが必要 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 最適化法のチューニング (1) 負のモル数 得られた組成が負になることがある 解くべき連立方程式が非正則になりやすい 行要素のほとんどがゼロ 計算の初期組成を平衡解に十分近付ける必要あり 局所的な極値へ収束する可能性 極少量の化学種の除去 化学ポテンシャルが発散する %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 最適化法のチューニング (2) 収束条件の設定 n_k+1 と n_k をどこまで接近させるか 低温で物性値が発散しないようにする 飽和蒸気圧を使う場合には注意 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 負のモル数 方程式系に n > 0 の条件は入ってない 何らかの条件を入れて組成が負になるのを防ぐ必要がある %% %filter "platex2eps.sh fu" \color[named]{Red} {\boldmath \begin{eqnarray} n_{k+1}^{\phi *} = n^{\phi}_{k} + \chi \{n^{\phi}_{k+1} - n^{\phi}_{k}\} \nonumber \end{eqnarray} } %endfilter %% %filter "platex2eps.sh fu-2" \color[named]{Red} {\boldmath \begin{eqnarray} \chi_{i} \equiv \left\{ \begin{array}{ll} \frac{n^{\phi}_{k}}{n^{\phi}_{k} - n^{\phi}_{k+1}} \times 0.9 & (n_{i}^{\phi} \leq 0) \\ 1 & (n_{i}^{\phi} > 0) \end{array} \right. \nonumber \end{eqnarray} } %endfilter %% %filter "platex2eps.sh nk" {\boldmath $n^{\phi *}_{k+1}$ } %endfilter %% 化学種のモル数が負になった場合 探索ベクトルの長さを短くする χはχ_i の最小値とする %cont %image "nk.eps" 0 250 250 1 %cont をより真に近い組成とみなす %center %image "fu.eps" 0 250 250 1 %image "fu-2.eps" 0 250 250 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 非正則/局所的な極値の回避 当該温度圧力において存在しうる化学種のみで計算 行要素がほとんどゼロにならないように 組成の 2 次方程式にならないケースを回避 1 成分 2 相系 2 相に含まれる化学種が全く同じ 平衡組成に十分近い組成を初期値とする %mark %cont, image "images/local-minimum.eps" 0 160 160 1 %size 3 %again %right 初期組成が平衡組成に十分に 近くないと, 最小値でなく 局所的な極値に収束する可能性有り %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 初期値の作り方 %size 2 平衡組成に近い組成を作成 1)系に含まれる化学種からある凝縮相を 1 つ選択. 気相とその凝縮相の 2 相からなる系を作る 2)前の圧力ステップでの平衡組成を凝縮相の初期値とする 気相の初期値は元素量から凝縮相の元素を除いたものから作成 もしも前の圧力ステップにおいてその凝縮相が存在しない場合は, 構成元素の適当な割合(3 割)を割り当てる 前のステップでの平衡組成が元素量よりも多い場合, 凝縮相内のモル比が保たれるように凝縮相のモル数を減らす 3) 上記の 1) の 2 相系の平衡状態を計算 当該温度・圧力において, 気相と平衡する凝縮相のモル数を計算 4) 上記の 1) -- 3) の操作を系に含まれる全ての凝縮相に対して実行 5) 上記の 4) から得られた平衡組成をそれぞれの凝縮相の初期組成と見なす. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 極少量の化学種の除去 化学ポテンシャルを発散させない 化学種が少なくなると log の項が発散 %% %filter "platex2eps.sh mu" \color[named]{Red} {\boldmath \begin{eqnarray} \mu_{i}^{\phi} = {\mu_{i}^{\phi}}^{\circ} + RT \ln{x_{i}^{\phi}} + \delta_{1\phi} RT \ln{\frac{p}{p_{0}}} \nonumber \end{eqnarray} } %endfilter %% %center %image "mu.eps" 0 250 250 1 %left 化学種の削除条件を設定 化学種に含まれる各々の元素比がマシンイプシロン程度になったら削除 モル比で設定すると, 1 成分で 1 相を形成する化学種を扱えない %% %filter "platex2eps.sh cut-mol" \color[named]{Red} {\boldmath \begin{eqnarray} \left(\frac{a_{ie} n_{i}}{\sum_{i} a_{ie} n_{i}} \right) \leq \varepsilon \nonumber \end{eqnarray} } %endfilter %% %center %image "cut-mol.eps" 0 250 250 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" 収束判定条件 収束 = 探索ベクトルが十分小さい 反復計算前後のモル数の相対誤差がマシンイプシロン程度 %% %filter "platex2eps.sh conv" \color[named]{Red} {\boldmath \begin{eqnarray} \left| \frac{{n_{i}^{\phi}}_{k+1} - {n_{i}^{\phi}}_{k}}{{n_{i}^{\phi}}_{k}} \right| \leq \varepsilon \nonumber \end{eqnarray} } %endfilter %% %center %image "conv.eps" 0 250 250 1 %left ギブス自由エネルギーの相対誤差で判定すると, 微小成分の収束が甘くなる 多量にある化学種の寄与だけで収束が判定されてしまう %% %filter "platex2eps.sh conv-2" \color[named]{Red} {\boldmath \begin{eqnarray} G = \sum_{i} {n_{i}^{\phi}}_{k} {\mu_{i}^{\phi}}_{k} \nonumber \end{eqnarray} } %endfilter %% %center %image "conv-2.eps" 0 250 250 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" エントロピーの与え方 エントロピーが負にならないようにする 飽和蒸気圧式(Antoine 式)をそのまま低温へ外挿するとマイナス無限大に発散 %% %filter "platex2eps.sh antoine" \color[named]{Red} {\boldmath \begin{eqnarray} \ln\{p^{*}(T)\} = A - \frac{B}{C + T} \nonumber \end{eqnarray} } %endfilter %% %center %image "antoine.eps" 0 250 250 1 %left エントロピーが 0 K で値が 0 となるようにする %image "images/NH3_off.gif" 0 100 100 1 %cont %image "images/NH3_on.gif" 0 100 100 1 %size 3 %center (左)蒸気圧式をそのまま利用, (右)飽和蒸気圧が 10 Pa 以下の領域は線形補間 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %nodefault %fore "red", size 8, back "darkblue", font "thick", vgap 0, ccolor "gray" %bgrad 0 0 128 0 1 "white" "white" "blue" "white" "white" "white" "white" "white" %center, fore "yellow", font "thick" 計算結果 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" Settings of calculation Pressure range %filter "platex2eps.sh region" \color[named]{Red} {\bf \boldmath \begin{tabular}{|l|r|r|r|r|}\hline & Jupiter & Saturn & Uranus & Neptune \\ \hline range of pressure (bar) & 0.1 -- 20 & 0.5 -- 30 & 0.1 -- 250 & 0.1 -- 600\\ \hline \end{tabular} } %endfilter %cont %image "region.eps" 0 200 200 1 Reference pressure and temperature for each planet %filter "platex2eps.sh ref" \color[named]{Red} {\bf \boldmath \begin{tabular}{|l|r|r|r|r|}\hline & Jupiter & Saturn & Uranus & Neptune \\ \hline pressure (bar) & 0.6 & 1.2 & 2.3 & 1.2 \\ temperature (K) & 150 & 143 & 101 & 79 \\ \hline \end{tabular} } %endfilter %cont %image "ref.eps" 0 200 200 1 %size 3 Adiabat for each planet is determined as the curve which pass the reference pressure and temperature given above These values are arbitrary given from observations %left Chemical species included in the calculation %filter "platex2eps.sh sosei" \color[named]{Red} {\bf \boldmath \begin{tabular}{|l|}\hline chemical species \\ \hline H$_2$(g), He(g), Ne(g), H$_2$O(g), CH$_4$(g), NH$_3$(g), H$_{2}$S(g) \\ H2O$_2$(l), NH$_3$(l), H2S(l), CH4(l), \\ H$_2$O(s), NH$_3$(s), H2S(s), CH4(s), NH$_{4}$SH(s) \\ \hline \end{tabular} } %endfilter %size 3 %cont %image "sosei.eps" 0 200 200 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" Settings of calculation (cont) Elemntal abundance H, He, Ne: 1 x solar C, N, O, S: %cont, fore "red" 1 x solar, 5 x solar, 10 x solar %fore "black" Three sets of abundances are used, since condensible volatile abundances are not well known %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" Adiabatic lapse rate of Jupiter. %center %image "ps/jupiter_dtdz.png" 0 100 100 1 %%%cont %% %%%cont %%%image "ps/H2.gif" 0 50 50 1 %left %size 4 Adiabatic lapse rate of Jupiter. red: 1 x solar, blue: 5 x solar, violet: 10 x solar %size 2 %size 3 Adiabatic lapse rate is about 2 K/km The degree of lapse rate with the depth is mainly caused by the increase of H2(g) specific heat %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" Static stability of Jupiter %center %image "ps/jupiter_stab.png" 0 100 100 1 %left %size 4 Static stability of Jupiter red: 1 x solar, blue: 5 x solar, violet: 10 x solar %size 2 %size 3 1x solar case, max value of static stability is 3 x 10^-5 10x solar case, max value of static stability is 7 x 10^-5. this value is close to that of the Earth %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" Static stability of Saturn %center %image "ps/jupiter_stab.png" 0 100 100 1 %cont %cont %image "ps/saturn_stab.png" 0 100 100 1 %left %size 4 static stability of Saturn and Jupiter red: 1 x solar, blue: 5 x solar, violet: 10 x solar %size 2 %size 3 if 10 x solar, the ststic stability of the Saturn is smaller than that of the Jupiter The difference of N^2 between Saturn and Jupiter is caused by magnitude of gravity %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" Static stability of Uranus %center %image "ps/jupiter_stab.png" 0 100 100 1 %cont %cont %image "ps/uranus_stab.png" 0 100 100 1 %left %size 4 static stability of Uranus and Jupiter red: 1 x solar, blue: 5 x solar, violet: 10 x solar %size 2 %size 3 the static stability by CH4(s) is the most strong, but the magnitude of it for 10 x solar is close to that of the Jupiter for 1 x solar %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" cf: Static stability and cloud density of Uranus %center %image "ps/uranus_stab.png" 0 100 100 1 %cont %cont %image "ps/uranus_cldns.png" 0 100 100 1 %left %size 4 static stability and cloud density of Uranus red: 1 x solar, blue: 5 x solar, violet: 10 x solar %size 2 %size 3 N^2 of CH4 condensation layer is larger than that of H2O, since H2O condensation layer is located in the deeper region and H2O mixing ratio per unit mass dry air (g/kg) is small %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" Conclusion We estimate static stability of Jovian planet atmospheres Calculation method We newly developed a calculation method to obtain atmospheric equilibrium composition by minimizing Gibbs free energy for a given arbitrary elemental abundance %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %page %bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" Conclusion (cont) Static stability N^2 of Jupiter in the water condensation layer ranges from 10^-5 for 1 x solar to 10^-4 for 10 x solar; the latter is close to that of the Earth N^2 of Jupiter's atmosphere is about 10 times as large as those of the other Jovian planet atmospheres because of its large gravity For Uranus and Neptune, N^2 of CH4 condensation layer is larger than that of H2O, since H2O condensation layer is located in the deeper region and H2O mixing ratio per unit mass dry air (g/kg) is small %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%page %%%bgrad 0 0 128 0 1 "blue" "white" "white" "white" "white" "white" "white" %%Remarks %% If we consider the effect of mean molecular weight gradient, static stability will increase by a factor of 2 to 3 %% %%%image "ps/seisou.gif" 0 100 100 1