スポンサーサイト

上記の広告は1ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

Scilabで数値積分:地球深部の密度と圧力

Scilabで数値積分: 固体の比熱ではintegrate関数を用いてデバイモデルから結晶固体の比熱を計算しました。
またScilabで乱数による関数の積分ではモンテカルロ法を用いて関数の積分を行いました。

これらはいずれも既知の関数f(x)について具体的な積分範囲、たとえばa≦x≦bでの積分値

\int_a^b f(x)\mathrm{d}x

を求める手法でした。

これらの方法では関数f(x)の形が式として分かっている必要があります。
しかしながら、離散的な数値データに対する積分値が計算したい場合もあります。

001_201311020335151bb.png

Fig.1: 地球の内部構造(世界初!地球中心部の超高圧高温状態を実現 ~ようやく手が届いた地球コア~より)


今回は地球内部の密度を表す数値表から、数値積分を使って重力加速度と圧力の深さ依存性を計算します。



標準地球モデル:PREM


地球内部の構造は、地震波観測からfig.1のような層構造をしていることが知られています。
もっとも有名な地球内部構造のモデルはDziewonski and Anderson (1981)によるPREM(Preliminary Reference Earth Model)と呼ばれるものです。

002_20131102033515160.png

Fig.2: 地球内部の地震波速度分布と密度分布。横軸は地球の中心からの距離(半径)。


PREMにはP波速度、S波速度、密度の深さごとの値、及びこれら3つの数値から計算できる重力加速度、圧力などが表として掲載されています。(参考:PREMの表地球を貫通するトンネル内の落下運動)

今回はPREMの密度から重力加速度と圧力を計算し、PREMの論文の中の値と比較して計算結果のチェックを行います。

重力加速度の計算


高校物理では重力加速度は約9.8[m/s]として計算をすると思います。
重力加速度はは大きな質量を持つ地球から万有引力によって引っ張られる力に起因しているため、地球内部では重力加速度の値が地球表面の値とは異なります。

地球の中心からr[m]での密度をρ(r)[g/m3]とするとr[m]より深い部分の質量M(r)[g]は

M(r) = 4 \pi \int_0^r u^2 \rho(u) \mathrm{d}u

重力加速度は、万有引力の法則より

g(r) = G\frac{M(r)}{r} = \frac{4\pi G}{r}\int_{0}^{r}u^2\rho(u)\mathrm{d}u

ここでGは万有引力定数G=6.67259×10-14[m3/s2/g]です。
PREMの数値データはPREM.txtに保存しfscanfMat等を用いてScilabから読み出しを行います。離散データの数値積分にはinttrapを用いました。

重力加速度計算の積分のところだけ抜き出すと、以下のようになります。

// *** 重力加速度の計算 ***
Gravity_num = zeros(Radius);
for i = 2:length(Gravity_num)
Gravity_num(i) = 4 * %pi * grav / (Radius(i) ^ 2) * inttrap(Radius(1:i), Radius(1:i) .^ 2 .* RHO(1:i));
end


圧力の計算


地球内部では非常に長い地質学的時間のなかで静水圧平衡が成り立っていると考えることが出来るため、圧力の計算が出来ます。

\frac{\mathrm{d}P}{\mathrm{d}r} = - \rho(r)g(r)

これを素直に積分すると以下の式が得られます。

\int\mathrm{d}P = - \int\rho(r)g(r)\mathrm{d}r

文字を置き換えて

\int\mathrm{d}Q = - \int\rho(u)g(u)\mathrm{d}u

地球表面rsでの圧力Ps(つまり大気圧)から地球深部rの圧力Pまで積分すると

\int_{P_s}^{P}\mathrm{d}Q = - \int_{r_s}^r\rho(u)g(u)\mathrm{d}u

大気圧はPs=1013.25[hPa]ですが、これは地球内部の圧力である数[GPa](数万気圧)と比べて充分小さいので左辺はほぼ地球深部の圧力Pと考えることが出来ます。

\int_{P_s}^{P}\mathrm{d}Q = P - P_s \simeq P

したがって圧力は以下のような積分で計算することが出来ます。

P = \int_{r}^{r_s}\rho(u)g(u)\mathrm{d}u

この積分を行う部分をScilabのプログラムから抜粋すると以下のようになります。

// *** 圧力の計算 ***
Pressure_num = zeros(Radius);
for i = 1:length(Pressure_num)-1
Pressure_num(i) = inttrap(Radius(i:$), RHO(i:$) .* Gravity_num(i:$));
end


計算結果


重力加速度と圧力の計算結果をFig.3に示します。

003_20131102033515738.png

Fig.3: 地球内部の重力加速度(上)と圧力(下)


赤の実線で示されたのがPREMの表にある値で、青の破線で示したのが今回積分から計算された値です。多少の誤差が発生しているようですが、おおよそ上手く計算できています。

全ソースコードはDeepEarth_sce.txtです。

clear;

// *** 物理定数 ***
grav = 6.67259e-14; // 万有引力定数 (m^3/s^2 g)

// *** PREMのテーブルを読み出し ***
X = fscanfMat('PREM.txt');
Radius = 1E3 * X(:,1); // 半径 (m)
Vp = X(:,2); // P波速度 (m/s)
Vs = X(:,3); // S波速度 (m/s)
RHO = 1E3 * X(:,4); // 密度 (g/m^3)
Ks = 1E12 * X(:,5); // 断熱体積弾性率(Pa)
Mu = 1E12 * X(:,6); // 剛性率(Pa)
Nu = 1E12 * X(:,7); //
Pressure = 1E12 * X(:,8); // 圧力(Pa)
Gravity = X(:,9); // 重力加速度 (m/s^2)

// *** 重力加速度の計算 ***
Gravity_num = zeros(Radius);
for i = 2:length(Gravity_num)
Gravity_num(i) = 4 * %pi * grav / (Radius(i) ^ 2) * inttrap(Radius(1:i), Radius(1:i) .^ 2 .* RHO(1:i));
end

// *** 圧力の計算 ***
Pressure_num = zeros(Radius);
for i = 1:length(Pressure_num)-1
Pressure_num(i) = inttrap(Radius(i:$), RHO(i:$) .* Gravity_num(i:$));
end

// *** 地震波速度と密度のプロット ***
scf(0)
// 地震波速度
xsetech([0,0,1,0.45]);
plot(1E-3 * Radius, 1E-3 * Vp,'-r');
plot(1E-3 * Radius, 1E-3 * Vs,'-b');
legend(["Vp","Vs"],1);
ylabel("Velocity (km/s)");
xlabel("Radius (km)");
// 密度
xsetech([0,0.45,1,0.45]);
plot(1E-3 * Radius, 1E-3 * RHO,'-r');
ylabel("Density (kg/s^3)");
xlabel("Radius (km)");

// *** 重力加速度と圧力のプロット ***
scf(1)
// 重力加速度
xsetech([0,0,1,0.45]);
plot(1E-3 * Radius,Gravity,'-r');
plot(1E-3 * Radius,Gravity_num,'--b');
legend(["Table","Numerical"],4);
ylabel("Acceleration of gravity (m/s^2)");
xlabel("Radius (km)");
// 圧力
xsetech([0,0.45,1,0.45]);
plot(1E-3 * Radius,1E-12 * Pressure,'-r');
plot(1E-3 * Radius,1E-12 * Pressure_num,'--b');
legend(["Table","Numerical"],1);
ylabel("Pressure (GPa)");
xlabel("Radius (km)");


台形積分とスプライン積分


今回使用したinttrapは、データ店の間を直線で補完する(要するに台形の面積を計算する)台形積分です。
もう少し高度な計算をするスプライン積分もintsplinを使って計算することが出来ます。
入力するパラメータも同じなので、多くの場合そのまま置き換えることが出来るでしょう。

ただし今回のPREMのようにデータが不連続となる点で複数の値を持っているような表の積分を行うときにはintsplinでは下記のようなエラーが出ます。

!--error 999 
splin: 入力引数 #1 の値が間違っています: (厳密に)昇順ではありません, もしくは +-inf を検出しました.


このような場合は、今回のように台形積分inttrapを使います。

関連エントリ




参考URL




付録


このエントリで使用したScilabのシミュレーション用ファイルを添付します。ファイル名末尾の".txt"を削除して、"_"を"."に変更すれば使えるはずです。(参考:ねがてぃぶろぐの付録)


参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。
スポンサーサイト

tag: Scilab 数値積分  マントル 

ダイヤモンド号で行く地底旅行

電子工作ブログであるねがてぃぶろぐの読者の方々には、興味のある方は少ないかもしれませんが、地球深部の物理・化学について、一般の方にも比較的分かり易く書かれている入舩徹男著ダイヤモンド号で行く地底旅行を紹介します。




フロンティアと探査機


探査機はやぶさが、小惑星イトカワから、その表面のサンプルを持ち帰った事は記憶に新しいと思います。こういったニュースを聞くと、人類の科学技術の進歩により、地球以外の惑星でさえも人類にとって、いずれ到達することの出来るフロンティアとなるだろうと感じます。

000_20111016172656.jpg
小惑星イトカワに着地する「はやぶさ」(想像図) (c) 池下 章裕
JAXAのウエブページより


一方で、人類未踏の領域があるのは、なにも地球の外だけではありません。
それは、地球の深部です。
宇宙開発と比べると、いささか地味に感じられるかもしれませんが、地球の内部は、宇宙以上に良く分かっていない場所です。

岡山大学理学部地球惑星科学科のウエブページには、高校生向けの研究紹介として地球の内部構造の図が掲載されています。



岡山大学の描いた地球の内部構造図では、地球内部で起きている色々な現象が表現されています。実は、小惑星イトカワに向かって探査機『はやぶさ』を送り出したのと同様に、地球内部を調べるための探査船『ちきゅう』を使って、海の底から穴を掘って調べようと言うプロジェクトが存在します。

002_20111016172656.jpg
地球深部探査船『ちきゅう』の写真。真ん中に建っている鉄塔からドリルを降下させて海の底をゴリゴリと掘り進んでいく。JAMSTECのウエブページより


陸ではなく、わざわざ海の底を掘るのは、一般的に、海洋地殻(岡山大の図で、水色で示したのが海、その下の茶色の部分)の方が大陸地殻(緑の部分)よりも薄く、上部マントル(薄黄色の部分)まで到達するために、掘らなければいけない量が少なくて済むからです。
しかしながら、2011年10月16日現在において、上部マントルまで穴を掘った人はいません。

ダイヤモンド号で行く地底旅行


したがって、地球内部の構造は、観測事実(地震波、磁場、地表で取れる岩石の化学組成、隕石の化学組成など)と実験や計算の比較を行うことから得られた推定です。

ダイヤモンド号で行く地底旅行は、こういった最新の研究を、あたかもダイヤモンドで出来た探査船に乗って見てきたかのように解説している本です。

岡山大学の地球断面図のコールドプリュームと書かれた部分に注目すると、(茶色で示した)海洋地殻物質が(薄い黄色の)上部マントルや(濃い黄色の)下部マントルの中まで沈み込んでいることが分かります。
『ダイヤモンド号』の旅は、この下降流(コールドプリューム)に乗って地球の深部まで進み、図の反対側に書いてある上昇流(ホットプリューム)に乗って帰ってくると言う旅程を進みます。

ダイヤモンド号で行く地底旅行は、専門書と大衆向けレーベルの中間ぐらいの難易度ですが、ダイヤモンド号からの『眺め』などは、とても詩的な表現がなされていて楽しい本です。

章構成は、下記の通りです。

  1. 未知への旅に出よう!
  2. 地球内部の運動と構造
  3. 沈み込むプレートに乗って
  4. 上部マントルの物質学
  5. 横たわるプレート
  6. 下部マントルをゆく
  7. 中心核へ
  8. 上昇するプリューム
  9. 帰還


地球科学カテゴリに関して


Twitterのプロフィールにもあるとおり私は電子工作が趣味の大学院生です。専攻は地球惑星科学と言うことになっています。

電子工作ブログに、地球惑星科学の話題もどうかと思いますが、最近は電子工作に割く時間も少なくなってきたし、科学・技術に興味がある方にとっては、地球科学も面白いはずだと信じてカテゴリ新設となりました。

実際のところ、高校地学は暗記教科の様相を呈していてあまり面白くありませんが、一方で、研究として行う地球科学は面白い学問だと思っています。
逆に、高校の物理・化学・生物学が面白いと感じた人は、生粋の物理・化学・生物学よりも地球科学のほうが専攻する学問として楽しめるのではないかとさえ感じます。

まあ、大抵の学問は何でも面白いのですが!

今回は、地球科学カテゴリのエントリの最初のものとして、私の研究分野で有名な愛媛大学の入舩徹男先生の本を紹介することにしました。

参考URL




参考文献/使用機器




フィードバック



にほんブログ村 その他趣味ブログ 電子工作へ

 ↑ 電子工作ブログランキング参加中です。1クリックお願いします。


コメント・トラックバックも歓迎です。 ↓      


 ↓ この記事が面白かった方は「拍手」をお願いします。

tag:  マントル ダイヤモンド 

FC2カウンター
カテゴリ
ユーザータグ

LTspiceAkaiKKRScilabmachikaneyamaKKRPSoCCPAOPアンプPIC強磁性モンテカルロ解析常微分方程式トランジスタodeインターフェース状態密度DOSecalj定電流PDS5022スイッチング回路半導体シェルスクリプト乱数レベルシフトHP6632A温度解析ブレッドボードI2CR6452A分散関係トランジスタ技術可変抵抗確率論数値積分反強磁性セミナー非線形方程式ソルバ絶縁バンドギャップ熱設計偏微分方程式バンド構造GW近似カオス三端子レギュレータLEDフォトカプラシュミットトリガISO-I2CA/DコンバータLM358USBカレントミラーTL431マフィンティン半径PC817C数値微分アナログスイッチ発振回路サーボ直流動作点解析74HC40532ちゃんねる標準ロジックチョッパアンプLDAアセンブラFFTbzqltyイジング模型ブラべ格子開発環境補間量子力学電子負荷BSchパラメトリック解析単振り子基本並進ベクトル熱伝導繰り返しGGAMaximaTLP621ewidthSMP相対論抵抗位相図ランダムウォークスピン軌道相互作用六方最密充填構造不規則合金FETコバルト失敗談QSGWcygwinスレーターポーリング曲線スイッチト・キャパシタラプラス方程式gfortranキュリー温度状態方程式条件分岐格子比熱TLP552LM555TLP521三角波NE555過渡解析FXA-7020ZRWriter509テスタ詰め回路MCUマントルダイヤモンドQNAPデータロガーガイガー管自動計測UPS井戸型ポテンシャルawk第一原理計算仮想結晶近似ブラウン運動差し込みグラフ平均場近似fsolve起電力熱力学OpenMPスーパーセル固有値問題最適化最小値VCAシュレディンガー方程式VESTAubuntu最大値面心立方構造PGAOPA2277L10構造非線型方程式ソルバ2SC1815fccフェルミ面等高線ジバニャン方程式ヒストグラム確率論マテリアルデザイン正規分布結晶磁気異方性interp1フィルタ初期値ウィグナーザイツ胞c/aルチル構造岩塩構造スワップ領域リジッドバンド模型edeltBaOウルツ鉱構造重積分SIC二相共存ZnOquantumESPRESSOCapSensegnuplotmultiplot全エネルギー固定スピンモーメントFSM合金ノコギリ波フォノンデバイ模型ハーフメタル半金属TeXifortTS-110不規則局所モーメントTS-112等価回路モデルパラメータ・モデルヒストグラムExcel円周率GimpトラックボールPC直流解析入出力文字列マンデルブロ集合キーボードフラクタル化学反応三次元Realforce縮退日本語最小二乗法関数フィッティング疎行列シンボル線種ナイキスト線図陰解法負帰還安定性熱拡散方程式EAGLECrank-Nicolson法連立一次方程式P-10クーロン散乱Ubuntu境界条件MBEHiLAPW軸ラベルトランスCK1026MAS830L凡例PIC16F785LMC662AACircuit両対数グラフ片対数グラフグラフの分割specx.f

最新コメント
リンク

にほんブログ村 その他趣味ブログ 電子工作へ
上記広告は1ヶ月以上更新のないブログに表示されています。新しい記事を書くことで広告を消せます。