HOME > natural science Laboratory > コンピュータ・シミュレーション講座 > 仮想物理実験室

VisualC++ と OpenGL を利用した仮想物理実験室
【2-8-2】ばね弾性力による単振動運動(1次元)のシミュレーション

文責:遠藤 理平 (2010年8月30日) カテゴリ:仮想物理実験室(268)

前節ばね弾性力による運動の運動方程式とアルゴリズムでは、 ばね弾性力がある場合の運動方程式の導出と、コンピュータでシミュレーションを行うためのアルゴリズムの導出を行ないました。 t=t_n[s]のときの加速度 a_n[m/s^2]は、

(2.8.1)
差分方程式

と、そのときの位置 x_n[m/s]によって決まる値となります。この加速度 a[m/s^2]から速度 v[m/s]、位置 x[m]を順番に計算することで、ばね弾性力のある場合のシミュレーションを行うことができます。前回のアルゴリズムが x軸方向の変位のみを仮定しているため、今回は初期条件として、初期位置(x_0, 0, 0), 初速度(0, 0, 0) と指定します。

今回は、上図のようにばねを縮めるように初期のボールを配置します。すると、ばね弾性力でボールは箱とは反対側に押し出されます。ある程度の伸びると今度は逆に、ばねが縮もうとして箱側に引っ張られます。そして、また最初の位置に戻ってきます。 その結果として、行ったり来たりの運動を繰り返すのです。 この運動のことは、単振動運動と呼ばれ、物理学では非常に重要な運動のひとつとなります。

ばね弾性力による運動のための計算アルゴリズム

動かない箱の位置とボールの位置・速度・加速度の初期値

動かない箱の一辺の長さ box_L 、 位置を (box_x, box_y, box_z) で指定します。 箱の位置は、箱の中心の座標なので、床にめり込まないように box_z = box_L/2.0 とします。 次にばねについてですが、ばねの自然長を L 、ばね定数を k とします。

k = 10.0;     //ばね定数
L = 30.0;     //ばねの自然長 
box_L = 10.0; //箱の一辺の長さ
box_x = -20.0;
box_y = 0.0;
box_z = box_L/2.0;  //箱が床にめり込まないように

ball_r = 4.0;//ボールの半径

// 位置・速度・加速度の初期値 
x = 0.0;
y = 0.0;
z = ball_r;
vx = 0.0;
vy = 0.0;
vz = 0.0;
ax = 0.0;
ay = 0.0;
az = 0.0;

速度・位置を逐次計算するアルゴリズム

//位置の計算
x  = x  + vx * dt;
y  = y  + vy * dt;
z  = z  + vz * dt;
//速度の計算
vx = vx + ax * dt;
vy = vy + ay * dt;
vz = vz + az * dt;
//加速度の計算
ax = - k * ((ball[i].x - box_x) - L);    //<-----ばね弾性力からの寄与
ay = 0.0;
az = 0.0;

このアルゴリズムでは誤差が大きく、時間と共に、単振動の振幅が大きくなってしまいます。 アルゴリズムの改善は、次に行ないます。

ばね弾性力による単振動のシミュレーション

つるつるの床に、動かない箱とばねで繋がったボールが置かれている状況のシミュレーションです。運動は、ばね弾性力のみによる単振動運動となります。

VisualC++ + OpenGL プログラミング

【0.2日目】仮想物理実験室の構築 (ver1.1)」をベースにプログラミングします。上記で解説した逐次計算を行うアルゴリズムをプログラミングします。

仮想物理実験室変数の定義

//////////////////////////////////////////////////////////////////////////
// 変数の定義
//////////////////////////////////////////////////////////////////////////
#define _BITMAP 0 //アニメーション作成用ビットマップの保存 0:しない 1:する
//--------------------------------------------------------
// 仮想物理実験室変数の定義
//--------------------------------------------------------
double t = 0.0;     //時刻
double dt= 0.01;    //時間刻み
int tn = 0;         //ステップ数
double k = 10.0;     //ばね定数
double L = 30.0;     //ばねの自然長 

//箱の場所
double box_L = 10.0; //箱の一辺の長さ
double box_x = -20.0;
double box_y = 0.0;
double box_z = box_L/2.0;  //箱が床にめり込まないように

// ボールの定義
double ball_r = 4.0;//ボールの半径
struct BALL {
  double x, y, z;
  double vx, vy, vz;
  double ax, ay, az;
};
const int N = 1;
BALL ball[N];
void SetUp(void){
  for(int i=0; i<N; i++ ){
    ball[i].x  = 0.0;
    ball[i].y  = 0.0;
    ball[i].z  = ball_r; //ボールがめり込まないように
    ball[i].vx = 0.0;
    ball[i].vy = 0.0;
    ball[i].vz = 0.0;
    ball[i].ax = 0.0;
    ball[i].ay = 0.0;
    ball[i].az = 0.0;
  }
}

ボールの位置「x,y,z」、速度「vx,vy,vz」、加速度「ax,ay,az」を 構造体「BALL」のメンバとします。「BALL」型の変数「ball1」を宣言と同時に初期値を設定しています。

速度・位置を逐次計算するアルゴリズム

//--------------------------------------------------------
// 計算と物体の描画
//--------------------------------------------------------
void Calculate(){
  t = dt * double(tn);
  //ofs1 << dt*tn << " " << ball1.x << " " << ball1.vx << " " << ball1.ax <

位置「ball1.x」「ball1.y」「ball1.z」
速度「ball1.vx」「ball1.vy」「ball1.vz」
加速度「ball1.ax」「ball1.ay」「ball1.az」
を逐次計算します。 先にも述べましたが、このアルゴリズムでは誤差が大きく、 時間と共に単振動の振幅が大きくなってしまいます。 アルゴリズムの改善は、次回行ないます。

VisualC++ と OpenGL を利用した仮想物理実験室

第0章 仮想物理実験室の構築

第1章 様々な運動

第2章 ニュートンの運動方程式

第3章 剛体の運動(エネルギー保存則と運動量保存則)

付録

  • 【A-1】参考文献
    ・(A-1-1)OpenGL について
    ・(A-1-2)VisualC++ について
    ・(A-1-3)物理シミュレーション
    ・(A-1-4)数値計算

未分類

力学

量子力学

波動論



▲このページのトップNPO法人 natural science トップ

関連記事

仮想物理実験室







▲このページのトップNPO法人 natural science トップ