大騒ぎ
妻の兄弟が勢揃いして遊びに来てます。
賑やかで楽しいもんです。
苦しい時の言葉
行く末に、宿をそことも定めねば、踏み迷うべき道もなきかな。
弾性波動方程式の差分法による数値シミュレーション
最近以下の論文を読んでいた。
A numerical experiment on wave propagation in an elastic quarter space, Yasuo Sato,
J. Phys. Earth,20,287-299,1972
かなり古い論文で、いまさらこの論文を読んでなにがしたいんだ？
と思われたそこのあなた！その通りです。
仕事の関係で地震関係をしているのだけれど、P波、S波といいながら、
見てもいないことを抜け抜けと言葉に出していることが自分の中で許せなかった。
そこで上記の論文をトレースしてみて、自分でシミュレーションをしてみようという訳です。
言語にはC言語を使いました、アニメーション作成にはgnuplotで時刻毎のPNGファイルをつなげてapng(アニメーションpng)にして、さらにアニメーションgifにするといいかにも回りくどいことをしました。
以下にプログラム(C言語)を示します。
間違い等あればぜひお知らせください。
#include "math.h"
//空間方法差分式の各係数
const double dtdxVp = sqrt(0.75);
const double alpminusbetadtdx = (1 - 0.55*0.55)*dtdxVp*dtdxVp;
//double dtdxVs = dtdxVp*dtdxVp*0.55*0.55; obviously bug 2017.8.1
const double dtdxVs = dtdxVp*0.55;
//u-空間方向の2階時間微分の右辺式
double uDifferenceOperator(int i, int j, double u[60][60], double w[60][60])
{
//空間方向の差分式の値をリターン
return dtdxVp * dtdxVp * (u[i + 1][j] - 2.0 * u[i][j] + u[i - 1][j])
+ alpminusbetadtdx*(1.0/4.0)*(w[i + 1][j + 1] + w[i - 1][j - 1] - w[i + 1][j - 1] - w[i - 1][j + 1])
+ dtdxVs*dtdxVs*(u[i][j + 1] - 2.0 * u[i][j] + u[i][j - 1]);
}
//w-空間方向の2階時間微分の右辺式
double wDifferenceOperator(int i, int j, double u[60][60], double w[60][60])
{
//空間方向の差分式の値をリターン
return dtdxVs * dtdxVs * (w[i + 1][j] - 2.0 * w[i][j] + w[i - 1][j])
+ alpminusbetadtdx*(1.0 / 4.0)*(u[i + 1][j + 1] + u[i - 1][j - 1] - u[i + 1][j - 1] - u[i - 1][j + 1])
+ dtdxVp*dtdxVp*(w[i][j + 1] - 2.0 * w[i][j] + w[i][j - 1]);
}
int main()
{
//x方向の変位(時刻 t)
double u0[60][60] = { 0.0 };
//z方向の変位(時刻 t)
double w0[60][60] = { 0.0 };
//x方向の変位(時刻 t+1)
double u1[60][60] = { 0.0 };
//z方向の変位(時刻 t+1)
double w1[60][60] = { 0.0 };
//x方向の変位(時刻 t+2)
double u2[60][60] = { 0.0 };
//z方向の変位(時刻 t+2)
double w2[60][60] = { 0.0 };
/***********************************************************/
/* 注意!!!
境界面x=0と境界面z=0は、配列ではi=1上、j=1上を実境界面とする。
すなわち、i=0やj=0は領域の外を意味するものとする。
*/
/***********************************************************/
//最大時刻
int TMAX = 70;
//時刻を進めて差分式を解く
for (int time = 1; time <= TMAX; time++)
{
///////////////////////////////////////////////////////////////////////////////////////////////
// 仮想境界(弾性体の外部)の値取得start
// これは真の境界条件を中心差分式を用いることで仮想境界に対する式を導出し、
// 仮想境界の値をメインの差分式に組み込む事で境界条件を間接的に満足するようにしている。
double calcCoefficient = 0.395; //詳しくはノートp4を参照。要はλ/(λ+2μ)の計算結果
//x=0上(仮想境界上)の値を計算
for (int j = 1; j < 59; j++)
{
w1[0][j] = w1[2][j] + u1[1][j + 1] - u1[1][j - 1];
u1[0][j] = u1[2][j] + calcCoefficient*(w1[1][j + 1] - w1[1][j - 1]);
}
//z=0上の値を計算
for (int i = 1; i < 59; i++)
{
//入力加振
if(time<=11)
{
if (i >= 11 && i <= 15)
{
w1[i][0] = w1[i][2] + calcCoefficient*(u1[i + 1][1] - u1[i - 1][1])
+ sin(2 * 3.141592*time / 12.0) * *1 / 2);
}
}
else
{
w1[i][0] = w1[i][2] + calcCoefficient*(u1[i + 1][1] - u1[i - 1][1]);
}
u1[i][0] = u1[i][2] + w1[i + 1][1] - w1[i - 1][1];
//x=0,z=0における条件 20170803 add start
u1[0][0] = u1[1][0];
w1[0][0] = w1[0][1];
//x=0,z=0における条件 20170803 add end
}
// 境界条件の指定end
///////////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////////
//メインの差分式
for(int i = 1; i < 59; i++)
{
for(int j = 1; j < 59; j++)
{
double udiff = uDifferenceOperator(i, j, u1, w1);
double wdiff = wDifferenceOperator(i, j, u1, w1);
u2[i][j] = 2 * u1[i][j] - u0[i][j] + udiff;
w2[i][j] = 2 * w1[i][j] - w0[i][j] + wdiff;
}
}
///////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////
//時刻毎のベクトル場をpngファイルに保存
char filename[32] = "";
sprintf(filename,"result%d.dat",time);
FILE *fp = fopen(filename, "w");
for (int c = 1; c < 50; c++)
{
for (int l = 1; l < 50; l++)
{
double x = c*0.3;
double y = -l*0.3;
fprintf(fp, "%f %f %f %f\n", x, y, u2[c][l], -w2[c][l]);
}
}
fclose(fp);
/////////////////////////////////////////////////////////////////
///////////////////////////////
//値更新
//u0 <- u1, w0 <- w1
for (int i = 0; i < 59; i++)
{
for (int j = 0; j < 59; j++)
{
u0[i][j] = u1[i][j];
w0[i][j] = w1[i][j];
}
}
//u1 <- u2, w1 <- w2
for (int i = 0; i < 59; i++)
{
for (int j = 0; j < 59; j++)
{
u1[i][j] = u2[i][j];
w1[i][j] = w2[i][j];
}
}
/////////////////////////////////
}
return 0;
}
*1:1 + cos(2 * 3.1415 * i / 12.0
