マントル対流数値シミュレーション概論 B 非デカルト座標系でのマントル対流の基礎方程式系 D 誤差関数

Appendix C 高温高圧下での物質の弾性的・熱的性質

(February 9, 2022)

以下の記述は [11] を参考にした。

C.1 Birch-Murnaghan の状態方程式

高圧の地球惑星内部では、物質は大きな静水圧縮を受ける。その場合にはもはや無限小歪の取り扱いでは用をなさず、有限歪の取り扱いが必要となる。 Birch-Murnaghan の状態方程式とは、Murnaghan の一般的な有限歪弾性論を援用して、Birch が導いたものである。以下では、オイラー流の記述を用いた導出を行う。

変形前 (圧力 $0$ ) の状態での任意の微小ベクトル $d\bm{X}$ が静水圧による等方的な変形を受けて $d\bm{x}$ になったとすると、これらの間には

d\bm{x}=(1+\alpha)d\bm{X}

(C.1)

なる関係がある。ここで $\alpha$ は定数で、圧縮のとき $\alpha<0$ である。これよりオイラー歪 (Almansiの歪テンソル) $\bm{A}$ は

d\bm{x}\cdot d\bm{x}^{\ast}-d\bm{X}\cdot d\bm{X}^{\ast}=\left[1-\dfrac{1}{(1+% \alpha)^{2}}\right]d\bm{x}\cdot d\bm{x}^{\ast}=d\bm{x}\cdot 2\bm{A}\cdot d\bm{% x}^{\ast}

となって、有限等方歪 $e$ は

e=\dfrac{1}{2}\left[1-\dfrac{1}{(1+\alpha)^{2}}\right]

(C.2)

と書ける。この定義では $e$ は圧縮のときに負となり、高圧下での圧縮を論ずる際には少々不便である。そこで $e$ と符号を反転させた $f\equiv-e$ を以下の議論では用いることにする。

まず、有限変形に伴なう密度変化を求めよう。変形前 (圧力0) の体積素片の体積を $dV_{0}$ 、密度を $\rho_{0}$ とし、圧力下でこれらが各々 $d V$ 、 $\rho$ になったとすると、変形の等方性(式(C.1)) から

\dfrac{\rho}{\rho_{0}}=\dfrac{dV_{0}}{dV}=\dfrac{1}{(1+\alpha)^{3}}=(1-2e)^{% \frac{3}{2}}=(1+2f)^{\frac{3}{2}}

(C.3)

あるいは

f=\dfrac{1}{2}\left[\left(\dfrac{\rho}{\rho_{0}}\right)^{\frac{2}{3}}-1\right]

(C.4)

なる関係を得る。

高圧下で有効な状態方程式を導くために、熱力学的な考察を加える。体積を $V$ 、圧力を $p$ 、温度を $T$ 、エントロピーを $S$ 、内部エネルギーを $U$ 、 Helmholtz の自由エネルギーを $F$ とすると、

p=-\left(\dfrac{\partial{U}}{\partial{V}}\right)_{S}=-\left(\dfrac{\partial{F}% }{\partial{V}}\right)_{T}

(C.5)

なる関係がある(A.1.4章を参照せよ)。静水圧下の圧縮では、 $U$ 及び $F$ はそれぞれ断熱的、等温的に蓄えられた歪エネルギーに相当すると考えられる。

歪エネルギー $W$ ( $U$ または $F$ )を歪 $f$ のベキ級数で表わそう。 $f=0$ で $W=0$ 及び $p=0$ であることを考えると、このベキ級数は

W=af^{2}+bf^{3}+\cdots

(C.6)

の如く、 $f$ の0次及び1次の項は含まれないことが分かる。式(C.6)を式(C.5)に代入し、式(C.3)を考慮すると、圧力 $p$ と歪 $f$ の関係は

p=-\dfrac{\partial{W}}{\partial{V}}=-\dfrac{\partial{W}}{\partial{f}}\dfrac{d{% f}}{d{V}}=\dfrac{1}{3V_{0}}(1+2f)^{\frac{5}{2}}(2af+3bf^{2}+\cdots)

(C.7)

と書ける。

式(C.7) に含まれる係数 $a$ 、 $b$ を、物理的に明確な意味を持つ量で書き換えてみよう。式(C.7)から、体積弾性率 (体積圧縮率の逆数、 A.1.7章を参照せよ) $K$ は

K=-V\dfrac{\partial{p}}{\partial{V}}=\dfrac{2}{9V_{0}}(1+2f)^{\frac{5}{2}}% \left[a+(7a+3b)f+\cdots\right]

(C.8)

と書ける。特に $f=0$ ( $p=0$ に相当) の場合を考えると、

K_{0}=\dfrac{2}{9V_{0}}a\quad\Rightarrow\quad a=\dfrac{9}{2}K_{0}V_{0}

(C.9)

を得る。ここで $K_{0}$ は $p=0$ での体積弾性率である。同様に、体積弾性率の圧力勾配 $K^{\prime}=dK/dp$ について、

\dfrac{\partial{K}}{\partial{f}}=K^{\prime}\dfrac{d{p}}{d{f}}

の両辺を計算すると、

\dfrac{2}{9V_{0}}(1+2f)^{\frac{3}{2}}[3(4a+b)+(49a+21b)f+\cdots]=K^{\prime}% \dfrac{1}{3V_{0}}(1+2f)^{\frac{3}{2}}[2a+(14a+6b)f+27bf^{2}+\cdots]

特に $f=0$ の場合には

\dfrac{2}{3V_{0}}(4a+b)=K^{\prime}_{0}\dfrac{2}{3V_{0}}a\quad\Rightarrow\quad% \dfrac{b}{a}=K^{\prime}_{0}-4

(C.10)

を得る。ここで $K^{\prime}_{0}$ は $p=0$ での体積弾性率の圧力微分値である。これらを式(C.7)に代入してやると、

	$\displaystyle p$	$\displaystyle=3K_{0}(1+2f)^{\frac{5}{2}}f\left[1+\dfrac{3}{2}(K^{\prime}_{0}-4% )f+\cdots\right]$
		$\displaystyle=\dfrac{3}{2}K_{0}\left[\left(\dfrac{\rho}{\rho_{0}}\right)^{% \frac{7}{3}}-\left(\dfrac{\rho}{\rho_{0}}\right)^{\frac{5}{3}}\right]\left\{1+% \dfrac{3}{4}(K^{\prime}_{0}-4)\left[\left(\dfrac{\rho}{\rho_{0}}\right)^{\frac% {2}{3}}-1\right]+\cdots\right\}$		(C.11)

を得る。式(C.11)を Birch-Murnagham の状態方程式と呼ぶ。 $K$ 及び $K^{\prime}$ に断熱条件での値を用いれば、断熱圧縮における圧力 $p$ と密度 $\rho$ の関係式を与える。また等温条件でのそれを用いれば、等温圧縮における $p$ と $\rho$ の関係式を与える。

C.2 Mie-Grüneisen の状態方程式

先に議論した Birch-Murnaghan の状態方程式では、温度が上がる効果は含まれていなかった。ここでは、温度が上がる効果を含めた状態方程式を検討する。

温度が上がるにつれて、物質中には格子の熱振動のモードがエネルギーの低いものから励起され、内部エネルギーが上がる。 $N$ 個の原子からなる系を考えれば、各原子は3方向の運動の自由度をもつから、系の自由度は $3N$ である²² 2 より厳密には系全体の自由度に相当する6を引くべきであろうが、 $N\gg 1$ から $3N-6\simeq 3N$ とした。。これは $3N$ 個の振動子 (フォノン) の集まりと考えられる。系の高温での性質を知るには、この $3N$ の振動子がどのようなエネルギー分布、つまり振動数の分布を持つかが問題となる。

量子力学及び統計力学によれば、振動数 $\nu$ の振動子のとり得るエネルギー準位は $nh\nu$ (ただし $n=0,1,2,\cdots$ ) であり、各準位をとる確率は $\exp(-nh\nu/kT)$ である。ここで $h$ はプランク定数、 $k$ はボルツマン定数である。従って、1つの振動子の平均エネルギー $u$ は

u=\dfrac{\displaystyle\sum_{n=0}^{\infty}nh\nu\exp\left(-\dfrac{nh\nu}{kT}% \right)}{\displaystyle\sum_{n=0}^{\infty}\exp\left(-\dfrac{nh\nu}{kT}\right)}=% \dfrac{h\nu\exp\left(-\dfrac{h\nu}{kT}\right)}{1-\exp\left(-\dfrac{h\nu}{kT}% \right)}

(C.12)

となる。ここで、第2式から第3式への変形には、分母が初項1、公比 $\exp(-\dfrac{h\nu}{kT})$ の無限等比級数であること、及び分子はそれを $1/kT$ で微分したものであることを用いた。系全体の内部エネルギー $U$ は、式(C.12)の $3N$ 個の振動子について和をとったものになるから、

U=U_{0}+\sum_{i=1}^{3N}\dfrac{h\nu_{i}\exp\left(-\dfrac{h\nu_{i}}{kT}\right)}{% 1-\exp\left(-\dfrac{h\nu_{i}}{kT}\right)}

(C.13)

ここで $i$ は振動子の番号、 $U_{0}$ は $T=0$ Kでの内部エネルギーである。 $U$ の具体的な表式を定めるには、 $3N$ 個の振動子の振動数 $\nu_{i}$ の分布を決める必要がある。

これを基にして状態方程式を導くために、熱力学的な考察を加える。内部エネルギー $U$ とHelmholzの自由エネルギー $F$ の間には

U=F+TS=F-T\left(\dfrac{\partial{F}}{\partial{T}}\right)_{V}=-T^{2}\left(\dfrac% {\partial(F/T)}{\partial{T}}\right)_{V}

なる関係がある(A章を参照のこと)。これより、

	$\displaystyle F$	$\displaystyle=-T\int\dfrac{U}{T^{2}}\text{d}T=-T\int\left[\dfrac{U_{0}}{T^{2}}% +\sum_{i=1}^{3N}\dfrac{h\nu_{i}\exp\left(-\dfrac{h\nu_{i}}{kT}\right)}{T^{2}% \left(1-\exp\left(-\dfrac{h\nu_{i}}{kT}\right)\right)}\right]\text{d}T$
		$\displaystyle=U_{0}+\sum_{i=1}^{3N}kT\ln\left[1-\exp\left(-\dfrac{h\nu_{i}}{kT% }\right)\right]$		(C.14)

の如く $F$ の表式を得る。式(C.5)に代入して圧力 $p$ を求めると、

	$\displaystyle p$	$\displaystyle=-\left(\dfrac{\partial{F}}{\partial{V}}\right)_{T}=-\dfrac{d{U_{% 0}}}{dV}-\sum_{i=1}^{3N}\dfrac{h\exp\left(-\dfrac{h\nu_{i}}{kT}\right)}{1-\exp% \left(-\dfrac{h\nu_{i}}{kT}\right)}\dfrac{d\nu_{i}}{dV}$
		$\displaystyle=-\dfrac{d{U_{0}}}{dV}-\dfrac{1}{V}\sum_{i=1}^{3N}\dfrac{h\nu_{i}% \exp\left(-\dfrac{h\nu_{i}}{kT}\right)}{1-\exp\left(-\dfrac{h\nu_{i}}{kT}% \right)}\dfrac{d\ln\nu_{i}}{d\ln{V}}$		(C.15)

を得る。ただし $U_{0}$ と $\nu_{i}$ は体積 $V$ だけの関数と仮定した。式(C.15)で、振動子の振動数が体積によって変化する割合を表わすパラメータ

\gamma_{i}=-\dfrac{{V}d\nu_{i}}{{\nu_{i}}dV}=-\dfrac{d\ln\nu_{i}}{d\ln{V}}

(C.16)

がどの振動子 $i$ についても $\gamma_{G}$ で一定であると仮定すると、式(C.15)は

\displaystyle p

\displaystyle=-\dfrac{d{U_{0}}}{dV}+\gamma_{G}\dfrac{U}{V}

(C.17)

と簡略化される。この式は Mie-Grüneisen の状態方程式と呼ばれる。ここで式(C.17)の右辺第1項 $-dU_{0}/dT=p_{0}$ は $T=0$ Kでの圧力を表わし、第2項 $\gamma_{G}{U}/{V}$ は温度上昇に起因する熱振動による圧力で、熱圧力とも呼ばれる。

なお、式(C.16)から定義される $\gamma_{G}$ は Grüneisen parameter と呼ばれる。これは振動子の振動数が体積変化によって変化する割合を表わしており、高温下における物質の性質を論じる際に重要な量である。式(C.17)を体積一定のもとに温度で微分すると、

\left(\dfrac{\partial{p}}{\partial{T}}\right)_{V}=\gamma_{G}\dfrac{1}{V}\left(% \dfrac{\partial{U}}{\partial{T}}\right)_{V}

これより

\gamma_{G}=\dfrac{V\left(\dfrac{\partial{p}}{\partial{T}}\right)_{V}}{\left(% \dfrac{\partial{U}}{\partial{T}}\right)_{V}}=\dfrac{V{\alpha}{K_{T}}}{C_{V}}=% \dfrac{V{\alpha}{K_{S}}}{C_{p}}

(C.18)

のように表わされる。ここで $\alpha$ は熱膨張率、 $K_{T}$ は等温体積弾性率、 $K_{S}$ は断熱体積弾性率、 $C_{V}$ は定積石比熱、 $C_{p}$ は定圧比熱であり、

\alpha=\dfrac{1}{V}\left(\dfrac{\partial{V}}{\partial{T}}\right)_{p}=-\dfrac{1% }{V}\left(\dfrac{\partial{V}}{\partial{p}}\right)_{T}\left(\dfrac{\partial{P}}% {\partial{T}}\right)_{V}=\dfrac{1}{K}\left(\dfrac{\partial{p}}{\partial{T}}% \right)_{V}

及び $K_{T}/K_{S}=C_{V}/C_{p}$ の関係式を用いた。式(C.18)は式(1.33)と等価であることに注意。即ち、 $\gamma_{G}$ はマクロな観測量で表わせていることになる。実際の鉱物の $K_{S}$ 、 $\alpha$ 、 $C_{p}$ 、 $V$ の測定データから $\gamma_{G}$ を求めると、その温度依存性は小さく、ある程度以上の高温ではほぼ一定とみなせる。また鉱物種の違いによる $\gamma_{G}$ の違いは小さく、およそ1から2の範囲に収まっている。

C.3 格子の熱振動に関するDebye モデル

振動子の振動数 $\nu_{i}$ がどのような分布をしているかを Debye モデルに従ってモデル化し、式(C.13)の内部エネルギー $U$ を見積もろう。 Debye モデルでは、 $N$ 個の原子からなる $3N$ 個の振動子を連続体の弾性波で置き換える。ただし、振動子の個数が $3N$ 個と有限であることは、弾性波の振動数 $\nu$ をある最大値 $\nu_{\text{max}}$ で打ち切ることで処理する。

波数ベクトル $\bm{k}=(k_{1},k_{2},k_{3})$ 、振動数 $\nu$ で伝わる平面波

A\exp[i(k_{1}x_{1}+k_{2}x_{2}+k_{3}x_{3}-2\pi\nu{t})]

を考える。この平面波が体積 $V=L^{3}$ の立方体の中を伝わるとしよう。周期境界条件のもとで存在しうる波は、波数 $(k_{1},k_{2},k_{3})$ が

0,\pm\dfrac{2\pi}{L},\pm\dfrac{4\pi}{L},\dfrac{6\pi}{L},\cdots

を満たすものに限られる。このことを $(k_{1},k_{2},k_{3})$ を座標とする波数空間で考えると、 $\dfrac{2\pi}{L}$ 間隔に分布する格子点の各々に対応する波が1つ存在することになる。弾性波の速度を $v$ として、振動数 $\nu$ 以下のモードの総数 $M$ は、半径 $k=\dfrac{2\pi\nu}{v}$ の体積に含まれるモードの総数であるから、

M=\left(\dfrac{L}{2\pi}\right)^{3}\dfrac{4}{3}\pi{k^{3}}=\left(\dfrac{L}{2\pi}% \right)^{3}\dfrac{4}{3}\pi\left(\dfrac{2\pi\nu}{v}\right)^{3}=\dfrac{4\pi{V}}{% 3}\dfrac{\nu^{3}}{v^{3}}

と書ける。ここで、単位体積あたりに存在しうる波の個数が $(\dfrac{L}{2\pi})^{3}$ であることを用いた。上式より、モードの数の分布関数 $D(\nu)$ は

\displaystyle D(\nu)=\dfrac{d{M}}{d{\nu}}=4\pi{V}\dfrac{nu^{2}}{v^{3}}

(C.19)

と与えられる。さらに弾性波には1方向につき1つの縦波( $v=v_{p}$ ) と2つの横波( $v=v_{s}$ )があることを用いて、式(C.19)の $v$ をこれらの平均値で置き換え、

\displaystyle D(\nu)=\dfrac{d{M}}{d{\nu}}=4\pi{V}{\nu^{2}}\left(\dfrac{1}{{v_{% p}}^{3}}+\dfrac{2}{{v_{s}}^{3}}\right)

(C.20)

としよう。ところで、全モード数が振動子の数 $3N$ となるようにするには、振動数 $\nu$ をある値 $\nu_{\text{max}}$ で打ち切らねばならない。式(C.20)を用いて積分すると

3N=\int_{0}^{\nu_{\text{max}}}D(\nu)d\nu=\dfrac{4}{3}\pi{V}\left(\dfrac{1}{{v_% {p}}^{3}}+\dfrac{2}{{v_{s}}^{3}}\right){\nu_{\text{max}}^{3}}

であるから、振動数の最大値 $\nu_{\text{max}}$ は

\displaystyle\nu_{\text{max}}=\left(\dfrac{9N}{4\pi{V}}\right)^{\frac{1}{3}}% \left(\dfrac{1}{{v_{p}}^{3}}+\dfrac{2}{{v_{s}}^{3}}\right)^{-\frac{1}{3}}

(C.21)

であり、式(C.20)と(C.21)より

\displaystyle D(\nu)=\dfrac{9N}{{\nu_{\text{max}}}^{3}}{\nu^{2}}

(C.22)

と与えられる。

モード数の分布が得られたので、式(C.13)の内部エネルギーを求めよう。式(C.12)と(C.22)より

\displaystyle U=\int_{0}^{\nu_{\text{max}}}{u}D(\nu)d\nu=\frac{9N}{{\nu_{\text% {max}}}^{3}}\int_{0}^{\nu_{\text{max}}}\dfrac{h\nu^{3}}{\exp\left(\dfrac{h\nu}% {kT}\right)-1}d\nu

(C.23)

であるが、ここで

\displaystyle\dfrac{h\nu}{kT}=\xi,\quad\dfrac{h\nu_{\text{max}}}{kT}=\dfrac{% \theta_{D}}{T}=x

(C.24)

とおけば、

\displaystyle U=9NkT\frac{1}{x^{3}}\int_{0}^{x}\dfrac{\xi^{3}}{e^{\xi}-1}d\xi=% 3NkT\,{f(x)}

(C.25)

と書ける。ここで

f(x)\equiv\frac{3}{x^{3}}\int_{0}^{x}\dfrac{\xi^{3}}{e^{\xi}-1}d\xi

とおいた。 $f(x)$ は高温 ( $1/x=T/\theta_{D}\to\infty$ ) で $f(x)\to 1$ となり、その結果 $U\to 3NkT$ である。

式(C.24)の $\theta_{D}$ は温度の次元をもつ量で、Debye 温度と呼ばれる。式(C.21)より

\displaystyle\theta_{D}=\dfrac{h\nu_{\text{max}}}{k}=\dfrac{h}{k}\left(\dfrac{% 9N}{4\pi{V}}\right)^{\frac{1}{3}}\left(\dfrac{1}{{v_{p}}^{3}}+\dfrac{2}{{v_{s}% }^{3}}\right)^{-\frac{1}{3}}

(C.26)

と書ける。 Debye 温度は、比熱や熱膨張などの熱的測定から得られるほか、弾性波速度 ( $v_{p}$ と $v_{s}$ ) からも得られる。造岩鉱物の Debye 温度はおよそ800K〜1000K程度の値であり、熱的に求めたものと弾性的に求めたものの両者はおおよそ一致する。

内部エネルギー $U$ の表式が式 (C.23) の通り得られたので、これを温度で微分してやると定積比熱 $C_{v}$ が求まる。

\displaystyle C_{v}=\left(\dfrac{\partial{U}}{\partial{T}}\right)_{V}=9Nk\frac% {1}{x^{3}}\int_{0}^{x}\dfrac{\xi^{4}{e^{\xi}}}{\left(e^{\xi}-1\right)^{2}}d\xi% =3Nk\,{c(x)}

(C.27)

ここで、

c(x)\equiv\frac{3}{x^{3}}\int_{0}^{x}\dfrac{\xi^{4}{e^{\xi}}}{\left(e^{\xi}-1% \right)^{2}}d\xi

とおいた。式(C.27) で $T\to\infty$ とすると $C_{v}\to 3Nk$ となり、古典的な Dulong-Petit の法則と一致する。