概要本研究では、粒子法の一種のMPS 法に基づいた、三次元の表面張力の新しいモデルを開発した。表面張力の数式((1)式)に基づいたモデルは、現在までに二次元の手法は開発されているが、それを三次元に拡張するのは難しい。そこで、三次元の計算モデルを開発した。また、その数値計算法をマグマオーシャンでの金属の分離の問題に適用した。

三次元の表面張力モデルの開発地球科学の中の流体系の多くは多相流体系である。多相流体では、物質の量の保存性を満たしつつ、相境界の移動を正確に計算機で計算することは難しい。この難点は差分法や有限要素法などのオイラー的な方法でなく、ラグランジュ的な方法である粒子法を用いることによって、解消される場合が多い。従って、本研究では粒子法の一種のMPS (Moving Particle Semi-implicit) 法を用いて、多相流体の計算コードを作成する。MPS 法は非圧縮流れに適した粒子法で、差分法などでは空間に固定されている計算点が、MPS 法では粒子に対応し、その計算点が流れによって移動する計算法である。

多相流体の中で、特にミクロな現象では界面に働く表面張力を正確に計算することが非常に大事である。表面張力は体積力として、次の式で与えられる。

(1)σ は表面張力係数(N/m)、n は界面の法線ベクトル、∇Σ は界面に沿った方向のgradient 演算子、δΣ は界面でのみ値を持つ表面デルタ関数である。ここで、∇Σ n は界面の平均曲率である。従って、表面張力の計算では、この三つの量、法線ベクトル、平均曲率、表面デルタ関数を計算する必要がある。

本研究では、(1) 式に基づいた三次元の表面張力モデルを開発した。現在までに存在するモデル(Nomura et al., (2001)[1]) は二次元にしか適用できず、三次元化も非常に困難であるので、新しいモデルを作成した。新しいモデルでは、法線ベクトルをカラーファンクション(θ = 0 at phase 0, θ = 1 at phase 1)の勾配から計算する。この勾配にはSPH 法の微分モデルを用いる。SPH 法は圧縮性流れに適した粒子法で、物理量を滑らかに記述できるのが強みである。ここでは、カラーファンクションのようなステップ状の関数を微分するときに安定に計算するために、SPH 法の微分モデルを用いた。また、曲率の計算にはMPS 法の微分モデルを用いた。SPH 法を用いると滑らかになりすぎ、大きな曲率の計算で大きな誤差が生じるためである。界面をシャープに記述するため表面デルタ関数には、Nomura らと同じMPS 法のモデルを用いた。この手法はアルゴリズムがややこしくなく、幾何学的なことを考えないですむので、三次元に適用しやすい。

計算例として、液滴の振動の現象を挙げる(図1)。正方形の液滴を静かに置くと、一定の周期で振動するという現象である。計算結果は理論値の周期(π1/4/√15)と一致している。

マグマオーシャンでの金属の液滴の計算本研究で開発したモデルの地球科学での応用例として、マグマオーシャンでの金属の分離の問題を取り扱った。金属とシリケイトがマグマオーシャン中で分離し、地球のマントルとコアを形成するというシナリオである。二次元と三次元の計算を行い、二次元の計算にはNomura らのモデルを、三次元の計算には本研究で開発したモデルを用いた。

このモデルは、1000km ほどの深さのマグマオーシャン中を、金属の液滴(-1cm)が、液体のシリケイトの中を密度差により降下して分離するというものである(Rubie et al., 2003[2])。なぜ、1cm になるかというと、表面張力と重力と抵抗力の三つが釣り合うサイズが1cm、ということである。液滴のサイズは化学物質のマントルとコアへの分配を考える上で非常に重要である。サイズが小さい、すなわち、体積に対し表面積が大きい場合、マグマオーシャンの下部の圧力下で、金属とシリケイトが化学平衡になっていたことになる。Rubie らによると、1cm のサイズだとすると、僅か200m の降下中に、化学平衡に達すると見積もられる。その場合、マグマオーシャン下部の圧力下で金属とシリケイトが化学平衡になるということである。

Rubie らの研究では、力の釣り合いを考えて、次元解析を行い結果を出したものである。従って、分裂や衝突の物理過程が入っていない。本研究はその物理過程をきちんと計算し、液滴のサイズや速度、サイズ分布を導出することが目的である。また、Rubie らの研究では安定な液滴のサイズを議論しただけだが、そのサイズに至るまでの過程を議論していない。例えば、初期状態として、非常に小さな液滴から分離が始まった場合や、既に形成されている大きな液滴から分離が始まった場合が考えられる。例えば、地球の形成物質と考えられるコンドライトに含まれるような小さい金属の液滴から分離が始まった場合、小さい液滴のストークス沈降速度は非常に遅いので、マグマオーシャンが固化するまでに、シリケイトから小さい液滴は分離しない、というようなことが考えられる。しかし、現在の上部マントルからはそのような金属の粒は見つかっていない。このように、数値シミュレーションをして実際にその現象を確かめる必要がある。

三次元の計算結果を一つ挙げる(図2)。白い粒は金属の計算粒子、周りの表示していない部分がシリケイトの計算粒子である。上下左右の境界は周期境界である。この計算例は、初期状態で液滴の数が一つの場合である。本計算により、液滴は合体と分裂を繰り返す非常に複雑な運動をすることがわかった。また、定常状態には10 秒程度で落ち着く。大きな液滴サイズから、計算を始めた場合も、小さなサイズから始めた場合も同じ定常状態に落ち着くことがわかった。定常状態での平均サイズや平均速度はRubie らの次元解析の値と比べて、オーダーは同じである。Rubie らの次元解析は平均値をよく再現しているということである。本研究により、次元解析で求めることの出来ない液滴のサイズ分布を算出できた。それによると、液滴のサイズに下限が存在することがわかった。これは、非常に小さい液滴が最初から存在していたとしても、合体により大きな液滴になり、すべての金属がシリケイトからマグマオーシャン中で速やかに分離したことを示唆する結果である。

コアとマントルの分離の過程で、どちらに熱エネルギーが多く分配されるか、という地球の熱史を考える上での重要な問題がある。本研究では、粘性散逸による温度の上昇も計算している。これは、上記の問題のマグマオーシャンでの分離過程に対する答えになる。計算結果によると、金属とシリケイトの温度差は1K にも満たない小さな温度差で抑えられることがわかった。これは無視できる値であり、マグマオーシャンでの分離過程では金属とシリケイトは、ほぼ等温で分離する考えられる。

計算では、粘性率をパラメータとして変化させた。マグマオーシャンでは粘性率が一番不確定だからである。粘性率を上昇させていくとき、粘性率が1.0 Pa・s 付近で、平均サイズの粘性に対する依存性が変化する(図3)。粘性率が1.0 Pa・s でのレイノルズ数を概算してみると、約3.4 であり、このあたりで慣性抵抗と粘性抵抗が逆転するようである。すなわち、粘性率が1.0 Pa・s よりも小さいところは、粘性がほとんど効かない慣性力が支配的な領域、粘性率が1.0 Pa・s よりも小さいところは、慣性がほとんど効かない粘性力が支配的な領域である。

図1: 液滴の振動のベンチマーク。(a) t = 0. (b) t = π1/4/2√15. (c) t = π1/4/√15. (d) t = 3π1/4/2√15.

図2: 三次元の金属の液滴の計算例:7.0 秒後(上下12cm、横6cm の領域、粘性率は1.0 Pa・s)

図3: 粘性率に対して、平均サイズをプロットした図。Morton 数は粘性率と表面張力係数、重力を含む無次元数。(二次元の計算結果)

本論文は、多相流体の振る舞いを理解する上で重要な表面張力を取り扱う新しいモデルを開発し、これを3次元多相流の数値計算に組み込むことに成功した。さらにこの新しいモデルと数値計算手法を用いて、初期地球におけるマグマオーシャン中での金属とシリケイトの分離過程のシミュレーションを行い、金属が沈降する様子を、速度、液滴のサイズ、温度・エネルギーなどの分布から詳細に検討した。その結果に基づいて、初期地球においてコアが形成される際の諸過程を議論し、その流れの様子、時間スケール、金属一シリケイト問の化学平衡・熱平衡の程度に制約を与えた。

第1章では、表面張力の重要性と本論文の概要が説明されている。特に3次元数値計算における表面張力の扱いについて、本研究でのアイデアと工夫、すなわちMoving-Particle-Selmi-implicit(MPS)methodとSmoothed Particle Hydrodynamics(SPH)を組み合わせて取り組む方針が提示されている。第2章では、表面張力を含む多相流れの支配方程式が整理、提示されている。特に、2相間の幾何学とそこに働く力を詳細に解析して、従来の2次元モデルの考え・手法を3次元に拡張することは困難であることを示し、3次元相境界を扱うための新たなパラメターの導入とそれを用いた定式化がなされた。これらの定式化に基づき、第3章において表面張力を含む3次元2相流の数値計算モデルが提示される。この新しいモデルでは、境界面の法線ベクトルを注目する物質の有無を表す関数の勾配から求めるが、その勾配をSPH法を用いて`滑らかに'求める一方、境界面の曲率はMPS法を用いてなるべくミクロな形状も考慮しながら求める。この工夫により、3次元の相境界を含む流れの計算が安定してかつ精度良く行えるようになった。本研究では、液滴の振動現象と渦輪(vortex ring)を再現することによってその精度を確認している(第3章、第3節)。

第4章では、本研究で開発された新しい数値計算モデルが、地球初期のマグマオーシャン中での金属とシリケイトの分離過程のシミュレーションに応用されている。同時に、広いパラメタースペースを探索するために、既存のモデルを用いて2次元シミュレーションも行われた。いずれの場合にも、初期の各相の量比、粘性、金属相の大きさ・分布を主要なパラメターとし、10cm×10cm程度の水平断面をもつ領域(主に液体シリケイトで占められる領域)中を金属の液滴がどのように落下するかが調べられた。その結果、2次元と3次元のいずれの場合にも、(1)落下する金属滴は合体と分裂を繰り返すこと、(2)さまざまな初期状態に関わらず速やかに(典型的には10秒程度で)定常的状態に達し、一定の速度分布、液滴サイズ分布を示すようになること、(3)定常的速度分布、液滴サイズ分布のパラメター依存性がレイノルズ数=1前後を境として変化することが明らかとなった。これまでに、次元解析を用いた比較的単純な議論から、マグマオーシャン中での金属液滴の代表的サイズが予測されており、今回の結果も、平均的サイズ(1cm程度)はこの予測に合うことが示された。同時に、これまで予測されていないサイズ分布について、下限(数mm程度)が存在し、比較的速やかに金属液滴が落下・分離することが明らかとなった。これらの落下過程でのエネルギー輸送・散逸も計算され、金属液滴とシリケイトの問の温度差は小さく(<1K)抑えられ、ほぼ等温で分離することも示された。第5章では、結果と結論がまとめられている。

本研究では、従来取り扱いが困難であった3次元の相境界とそこでの力学に正面から取り組み、これを多相流の数値計算モデルに組み込むことに成功した。さらに地球の進化を決める重要なプロセスの一つであるコア形成過程と、そこでの素過程であるマグマオーシャン中での金属とシリケイトの分離過程解明にこの新しいモデルを応用し、上記のような新たな知見を得た。本研究で開発された新しいモデルは、金属一シリケイトの分離過程のみならず、地球・惑星に普遍的に見出されるさまざまな多相流の関与する現象に応用可能であるという一般性を有する。

なお、本論文第3章は、S.Labrosseとの共同研究、第4章はS.Labrosse及び栗田敬との共同研究であるが、論文提出者が主体となってモデル開発、シミュレーション、解析を行ったもので、論文提出者の寄与が十分であると判断する。

したがって、博士(理学)の学位を授与できると認める。