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

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

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

(2-10-1)【2-10-1】ばね弾性力+重力による運動の運動方程式(3次元)とアルゴリズムでは、 ばねの変位が xyz平面内の場合で、ばね弾性力重力がボールに作用する場合の運動方程式の導出と、コンピュータでシミュレーションを行うためのアルゴリズムの導出を行ないました。 t=t_n[s]のときの加速度 a_n[m/s^2]は、

(2.10.1)
差分方程式

となり、そのときのボールの位置 r_n = (x_n, y_n, z_n) によって決まる値となります。 L_n は、原点からボールまでの距離、

(2.10.2)
差分方程式

です。コンピュータシミュレーションでは、 初期条件として、l_0 = 20.0[m]とし、10個のボールの質量 m [kg]を変化させて運動の違いを確認します。

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

支柱の位置とボールの位置・速度・加速度の初期値

ばね定数を k、自然長 L[m] のばねの片側を固定した動かない支柱を原点に置き、ばねの反対側をボールつけ、x-y 平面上に伸ばします。時刻 t=t_n[s] のときのボールの位置を r_n とします。

k = 10.0;    //ばね定数
L = 30.0;    //ばねの自然長
r = 20.0;    //ボールの初期半径
g =9.8
ball_r = 4.0;//ボールの半径

//支柱の定義
poll_x = 0.0;
poll_y = 0.0;
poll_z = 60.0;
poll_r = 5.0;  //半径
poll_h = 20.0; //高さ

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

//位置の計算
x  = x  + vx * dt;
y  = y  + vy * dt;
z  = z  + vz * dt;

l = sqrt(pow(ball[i].x,2)+pow(ball[i].y,2)+pow(poll_z-ball[i].z,2)); //<-----支柱からボールまでの距離

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

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

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

動かない支柱にばねで繋がったボールが運動するシミュレーションです。運動は、ばね弾性力と重力による運動となります。

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

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

ばねの定義

ばねの定義については、ばね弾性力による運動の運動方程式(1次元)とアルゴリズムをご覧下さい。

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

//////////////////////////////////////////////////////////////////////////
// 変数の定義
//////////////////////////////////////////////////////////////////////////
#define _BITMAP 0 //アニメーション作成用ビットマップの保存 0:しない 1:する
ofstream ofs1( "ball1.data" );
//--------------------------------------------------------
// 仮想物理実験室変数の定義
//--------------------------------------------------------
double t = 0.0;     //時刻
double dt= 0.05;    //時間刻み
int tn = 0;         //ステップ数
double k = 10.0;    //ばね定数
double L = 30.0;    //ばねの自然長
double l = 20.0;    //ばねの初期長
double g =9.8;      //重力加速度
double radian, radian_phi, radian_theta; 
double theta;

//支柱の定義
double poll_x = 0.0;
double poll_y = 0.0;
double poll_z = 60.0;
double poll_r = 5.0;  //半径
double poll_h = 20.0; //高さ

// ボールの定義
double ball_r = 4.0;//ボールの半径
struct BALL {
  double m;
  double x, y, z;
  double vx, vy, vz;
  double ax, ay, az;
};
const int N = 10;
BALL ball[N];
void SetUp(void){
  for(int i=0; i<N; i++ ){
    radian = 2.0 * PI * double(i)/double(N);
    ball[i].m  = 2.0*double(i+1)/double(N); //ボールの質量
    ball[i].x  = l * cos(radian);
    ball[i].y  = l * sin(radian);
    ball[i].z  = poll_z; //ボールがめり込まないように
    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 <<endl;
  for(int i=0; i<N; i++){
    //位置の算出
    ball[i].x = ball[i].x + ball[i].vx * dt;
    ball[i].y = ball[i].y + ball[i].vy * dt;
    ball[i].z = ball[i].z + ball[i].vz * dt;

    l = sqrt(pow(ball[i].x,2)+pow(ball[i].y,2)+pow(poll_z-ball[i].z,2)); // 支柱からボールまでの距離

    //加速度の算出
    ball[i].ax = - k/ball[i].m * (1.0 - L/l) * ball[i].x;
    ball[i].ay = - k/ball[i].m * (1.0 - L/l) * ball[i].y;
    ball[i].az = - k/ball[i].m * (1.0 - L/l) * (ball[i].z-poll_z) -g;
    //速度の算出
    ball[i].vx = ball[i].vx + ball[i].ax * dt;
    ball[i].vy = ball[i].vy + ball[i].ay * dt;
    ball[i].vz = ball[i].vz + ball[i].az * dt;
  }
 tn++;
}
void DrawStructure(){
  for(int i=0; i<N; i++){
    glPushMatrix();
      glMaterialfv(GL_FRONT, GL_AMBIENT, ms_silver.ambient);
      glMaterialfv(GL_FRONT, GL_DIFFUSE, ms_silver.diffuse);
      glMaterialfv(GL_FRONT, GL_SPECULAR, ms_silver.specular);
      glMaterialfv(GL_FRONT, GL_SHININESS, &ms_silver.shininess);
      glTranslated(ball[i].x, ball[i].y, ball[i].z); //平行移動値の設定
      glutSolidSphere(ball_r, 20, 20);            //引数:(半径, Z軸まわりの分割数, Z軸に沿った分割数)
    glPopMatrix();
  //バネ
    glPushMatrix();
      glMaterialfv(GL_FRONT, GL_AMBIENT, ms_yellow_rubber.ambient);
      glMaterialfv(GL_FRONT, GL_DIFFUSE, ms_yellow_rubber.diffuse);
      glMaterialfv(GL_FRONT, GL_SPECULAR, ms_yellow_rubber.specular);
      glMaterialfv(GL_FRONT, GL_SHININESS, &ms_yellow_rubber.shininess);  
      drowSolidSpring(box_x, box_y, box_z, (ball[i].x, (ball[i].y,(ball[i].z); //引数:(始点x, 始点y,  始点z, 終点x, 終点y, 終点z )
    glPopMatrix();
  }
  //軸(円錐)
  glPushMatrix();
      glMaterialfv(GL_FRONT, GL_AMBIENT, ms_copper.ambient);
      glMaterialfv(GL_FRONT, GL_DIFFUSE, ms_copper.diffuse);
      glMaterialfv(GL_FRONT, GL_SPECULAR, ms_copper.specular);
      glMaterialfv(GL_FRONT, GL_SHININESS, &ms_copper.shininess);
      glTranslated(poll_x, poll_y, poll_z);
      glutSolidCone(5.0,20.0,20,20);//引数:(半径, 高さ, Z軸まわりの分割数, Z軸に沿った分割数)
  glPopMatrix();
}

位置「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 トップ