062_Vibration1.mw

機械振動論 

1質点の振動 

振動論はあらゆる場面で重要な要素であり、この基礎理論を理解することはエンジニアにとって基本的なものです。 

ここでは、Mapleを使ってわかり易く振動論を勉強しようという試みにチャレンジしてみました。振動論の多くの教科書的な本が 

ありますが、Mapleを使えばインターネット時代の電子テキストとして有効な新しいツールを提供できます。 

> restart:
 

restartはMapleを初期化するコマンドです。これによりメモリー上に残っていたデータなどがクリヤーされます。 

> with(plots):with(plottools):
 

Warning, the name changecoords has been redefined 

Warning, the assigned name arrow now has a global binding 

これでプロット用のパッケージを呼び出しています。アニメーションなどで必要なので呼び出しておきます。 

アニメーション 

まずはアニメーションにより基本的な振動をアニメーション表示することにより視覚的に理解してみましょう。 

トップにあるアニメーションのように地面の上を一輪車が走っている状態を表現します。 

まず地面はplot3dで表現します。車輪はシリンダーで表し、これを回転(rotate)とシフト(translate)を使って地面の上に置きます。 

> p:=[]:
 

pはこれからアニメーションのフレームを追加していくための変数で、最初は空にしておきます。 

コロン(:=)は代入であり等号(=)と区別されます。 

> r:=1:
 

rは車輪の半径です。1にセットしておきます。 

> c:=cuboid([-0.2,-0.2,-0.2],[0.2,0.2,1.2],color=green):
 

cはタイヤにつながっている部分のオブジェクトを直方体(cuboid)で作成しておきます。 

対角の2点の座標で定義します。色は緑(green)にしておきます。 

これをあとでtranslateでシフト(平行移動)することで動かします。 

> M:=cuboid([-0.5,-0.5,0],[0.5,0.5,1],color=red):
 

Mは上に載っている物体を現す立方体です。色は赤にします。 

 

次に地面の形を現す関数g(t)を定義しましょう。 

> g:=t->sin(2*Pi*t)/10;
 

proc (t) options operator, arrow; 1/10*sin(2*Pi*t) end proc 

このように->を使って簡単に関数が定義できます。 

この関数をプロットしてみましょう。 

> plot(g(t),t=0..5);
 

Plot 

以下のforループで時間(tが0から1秒まで0.02秒刻み)ごとにアニメーションのフレームを作成します。 

スプリングはtubeplotを使って表現します。 

スプリングの長さをlで表します。 

hは車輪の中心の高さでHはその上の質点(赤)の中心の高さの変数です。 

> for t from 0 by 0.02 to 1 do
 h:=g(t):
 pp:=plot3d([x,y,g(t+x/4)],x=-Pi/2..Pi/2,y=-1..1,grid=[15,12],color=gray):
 pp:=pp,translate(rotate(cylinder([0,0,0],r,0.2,color=blue),Pi/2,t*2*Pi,0),0,-0.1,h+r):
 pp:=pp,translate(c,0,0,h+r):
 H:=sin(2*Pi*t)/4:
 l:=H-h+1:
 tt:=tubeplot([cos(u*8*Pi)/4,sin(u*8*Pi)/4,u*l],u=0..1,radius=0.05,color=yellow):
 pp:=pp,translate(tt,0,0,h+2*r+0.2):
 pp:=pp,translate(M,0,0,h+2*r+0.2+l):
 p:=p,display(pp):
end:
 

できたフレームをアニメーション表示します。 

> display(p,insequence=true,scaling=constrained);
 

Plot 

このモデルは単純化して1質点のモデルとしてシミュレーションすることができます。 

地面の変化が入力になりそれに対して上の質点の振動がその応答となります。 

次にこのモデルの運動方程式を考えてみましょう。 

> M:='M';g:='g';t:='t';c:='c';
 

M 

g 

t 

c 

これはさきほど定義した関数gおよび代入した変数を初期化し直しています。 

一旦代入されたものは式の中で変数として扱われませんので初期化する必要があります。 

>
 

数値シミュレーション ケース1 

数式をみても複雑で良くわかりませんので具体的な値を入れて数値シミュレーションを行いましょう。 

 

まず質量は1(Kg)としましょう。 

> M:=1;
 

1 

重力加速度はMKSでは9.8です。 

> G:=9.8;
 

9.8 

まずは簡単に、減衰定数cは0にします。 

> c:=0;
 

0 

バネ定数kは100とします。 

> k:=100;
 

100 

地面は平面で変化なし(平面)とします。 

> g:=t->0;
 

proc (t) options operator, arrow; 0 end proc 

これで運動方程式は次のようになります。 

> M*(diff(x(t),t$2)+G)+c*diff(x(t),t)+k*(x(t)-g(t))=0;
 

(diff(diff(x(t), t), t))+9.8+100*x(t) = 0 

これをdsolveで解いてみましょう。 

> dsolve({%,x(0)=0,D(x)(0)=0},x(t));
 

x(t) = -49/500+49/500*cos(10*t) 

x(0)=0は高さの初期値(t=0での値)を与えています。 

またD(x)はx(t)の一階微分でD(x)(0)=0は初速(t=0での速度の初期値)を0にしています。 

このように初期値を与えその後の変化を求める微分方程式の問題を初期値問題と呼びます。 

 

右辺を取り出しunapplyで関数にします。rhsは右辺(right hand side)を取り出すコマンドです。 

> ans:=unapply(rhs(%),t);
 

proc (t) options operator, arrow; -49/500+49/500*cos(10*t) end proc 

unapplyは式を関数化することができるたいへん便利な関数です。 

 

これをプロットしましょう。 

> plot(ans(t),t=0..3);
 

Plot 

バネ定数kの値を変えてやってみましょう。 

> k:=200;
 

200 

> M*(diff(x(t),t$2)+G)+c*diff(x(t),t)+k*(x(t)-g(t))=0;
 

(diff(diff(x(t), t), t))+9.8+200*x(t) = 0 

> dsolve({%,x(0)=0,D(x)(0)=0},x(t));
 

x(t) = -49/1000+49/1000*cos(10*2^(1/2)*t) 

右辺を取り出しunapplyで関数にします。 

> ans:=unapply(rhs(%),t);
 

proc (t) options operator, arrow; -49/1000+49/1000*cos(10*2^(1/2)*t) end proc 

これをプロットしましょう。 

> plot(ans(t),t=0..3);
 

Plot 

振動数が変わったのが解ります。 

このようにバネ定数は振動数に影響を与えます。同時に数値をみれば振幅も小さくなっていることが解ります。 

kの値が大きくなることはバネが硬くなりあまり振動しなくなると同時に振動数(周波数)が高くなります。 

>
 

数値シミュレーション ケース2 

今度はこれに減衰定数の効果を入れてみましょう。 

 

減衰定数cを0.1にしてみます。 

> c:=0.1;
 

.1 

これで運動方程式は次のようになります。 

> M*(diff(x(t),t$2)+G)+c*diff(x(t),t)+k*(x(t)-g(t))=0;
 

(diff(diff(x(t), t), t))+9.8+.1*(diff(x(t), t))+200*x(t) = 0 

これを解いてみましょう。 

> dsolve({%,x(0)=0,D(x)(0)=0},x(t));
 

x(t) = 49/79999000*exp(-1/20*t)*sin(1/20*79999^(1/2)*t)*79999^(1/2)+49/1000*exp(-1/20*t)*cos(1/20*79999^(1/2)*t)-49/1000
x(t) = 49/79999000*exp(-1/20*t)*sin(1/20*79999^(1/2)*t)*79999^(1/2)+49/1000*exp(-1/20*t)*cos(1/20*79999^(1/2)*t)-49/1000
 

右辺を関数にします。 

> ans:=unapply(rhs(%),t);
 

proc (t) options operator, arrow; 49/79999000*exp(-1/20*t)*sin(1/20*79999^(1/2)*t)*79999^(1/2)+49/1000*exp(-1/20*t)*cos(1/20*79999^(1/2)*t)-49/1000 end proc
proc (t) options operator, arrow; 49/79999000*exp(-1/20*t)*sin(1/20*79999^(1/2)*t)*79999^(1/2)+49/1000*exp(-1/20*t)*cos(1/20*79999^(1/2)*t)-49/1000 end proc
 

プロットしてみましょう。 

> plot(ans(t),t=0..5);
 

Plot 

減衰していることが確認できます。 

>
 

数値シミュレーション ケース3 

今度は地面の変化を入れてみましょう。 

> g:=t->sin(t)/5;
 

proc (t) options operator, arrow; 1/5*sin(t) end proc 

これで運動方程式は次のようになります。 

> M*(diff(x(t),t$2)+G)+c*diff(x(t),t)+k*(x(t)-g(t))=0;
 

(diff(diff(x(t), t), t))+9.8+.1*(diff(x(t), t))+200*x(t)-40*sin(t) = 0 

これを解いてみましょう。 

> dsolve({%,x(0)=0,D(x)(0)=0},x(t));
 

x(t) = -15725555051/316804119899000*exp(-1/20*t)*sin(1/20*79999^(1/2)*t)*79999^(1/2)+194444949/3960101000*exp(-1/20*t)*cos(1/20*79999^(1/2)*t)-49/1000-400/3960101*cos(t)+796000/3960101*sin(t)
x(t) = -15725555051/316804119899000*exp(-1/20*t)*sin(1/20*79999^(1/2)*t)*79999^(1/2)+194444949/3960101000*exp(-1/20*t)*cos(1/20*79999^(1/2)*t)-49/1000-400/3960101*cos(t)+796000/3960101*sin(t)
x(t) = -15725555051/316804119899000*exp(-1/20*t)*sin(1/20*79999^(1/2)*t)*79999^(1/2)+194444949/3960101000*exp(-1/20*t)*cos(1/20*79999^(1/2)*t)-49/1000-400/3960101*cos(t)+796000/3960101*sin(t)
 

右辺をunapplyにより関数にします。関数の名前はansです。 

> ans:=unapply(rhs(%),t);
 

proc (t) options operator, arrow; -15725555051/316804119899000*exp(-1/20*t)*sin(1/20*79999^(1/2)*t)*79999^(1/2)+194444949/3960101000*exp(-1/20*t)*cos(1/20*79999^(1/2)*t)-49/1000-400/3960101*cos(t)+796...
proc (t) options operator, arrow; -15725555051/316804119899000*exp(-1/20*t)*sin(1/20*79999^(1/2)*t)*79999^(1/2)+194444949/3960101000*exp(-1/20*t)*cos(1/20*79999^(1/2)*t)-49/1000-400/3960101*cos(t)+796...
proc (t) options operator, arrow; -15725555051/316804119899000*exp(-1/20*t)*sin(1/20*79999^(1/2)*t)*79999^(1/2)+194444949/3960101000*exp(-1/20*t)*cos(1/20*79999^(1/2)*t)-49/1000-400/3960101*cos(t)+796...
 

プロットしてみましょう。 

> plot({g(t),ans(t)},t=0..10);
 

Plot 

赤は地面の高さで緑が質点の高さです。 

この地面のように強制的に振動させようとすることを強制振動と呼びます。 

 

これをアニメーション表示しましょう。 

> p:=[]:
 

pはアニメーションのフレームを追加していくための変数で、最初は空にしておきます。 

> r:=1:
 

rは車輪の半径です。 

> c:=cuboid([-0.2,-0.2,-0.2],[0.2,0.2,1.2],color=green):
 

cはタイヤにつながっている部分を直方体(cuboid)で作成しておきます。 

対角の2点の座標で定義します。色は緑(green)にします。 

これをtranslateでシフトしてアニメーションします。 

> M:=cuboid([-0.5,-0.5,0],[0.5,0.5,1],color=red):
 

Mは上に載っている物体を現す立方体です。色は赤にします。 

 

以下のforループで時間(tが0から1秒まで)ごとにアニメーションのフレームを作成します。 

スプリングはtubeplotを使って表現します。スプリングの長さをlで表します。 

hは車輪の中心の高さでHは上の質点(赤)の中心の高さです。 

質点の高さはdsolveによる解のansを使用します。 

> for t from 0 by 0.05 to 5 do
 h:=g(t):
 pp:=plot3d([x,y,g(x+t)],x=-Pi/2..Pi/2,y=-1..1,grid=[15,12],color=gray):
 pp:=pp,translate(rotate(cylinder([0,0,0],r,0.2,color=blue),Pi/2,t*2*Pi,0),0,-0.1,h+r):
 pp:=pp,translate(c,0,0,h+r):
 H:=ans(t):
 l:=H-h+1:
 tt:=tubeplot([cos(u*8*Pi)/4,sin(u*8*Pi)/4,u*l],u=0..1,radius=0.05,color=yellow):
 pp:=pp,translate(tt,0,0,h+2*r+0.2):
 pp:=pp,translate(M,0,0,h+2*r+0.2+l):
 p:=p,display(pp):
end:
 

できたフレームをアニメーション表示します。 

> display(p,insequence=true,scaling=constrained);
 

Plot 

>
 

ステップ応答 

地面が一段の階段状(ステップ状)なっている場合のシミュレーションをしましょう。 

 

時間の変数tをリセットしておきます。 

> t:='t';
 

t 

> M:=1;c:=0.3;
 

1 

.3 

質量は1、粘性減衰定数は0.3とします。 

 

地面の変化を入れてみましょう。 

 

階段関数にはHeaviside関数を使用します。これは超関数(δ関数の積分)です。 

超関数がわからなくても以下のプロットのような階段上の関数だと思えば問題ありません。 

非常に小さい値epsを使用します。これは超関数Heavisideの階段の点での値が定義されないので 

わずかにずらすためです。 

> eps:=1e-10;
 

0.1e-9 

ステップ関数を定義します。 

> g:=t->0.1*Heaviside(t-0.5+eps);
 

proc (t) options operator, arrow; .1*Heaviside(t-.5+eps) end proc 

t=0.5(+eps)の点でステップが発生します。 

> plot(g(t),t=0..4);
 

Plot 

これで運動方程式は次のようになります。 

> M*(diff(x(t),t$2)+G)+c*diff(x(t),t)+k*(x(t)-g(t))=0;
 

(diff(diff(x(t), t), t))+9.8+.3*(diff(x(t), t))+200*x(t)-20.0*Heaviside(t-.4999999999) = 0
(diff(diff(x(t), t), t))+9.8+.3*(diff(x(t), t))+200*x(t)-20.0*Heaviside(t-.4999999999) = 0
 

これを解いてみましょう。 

> dsolve({%,x(0)=0,D(x)(0)=0},x(t));
 

x(t) = 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*79991^(1/2)*Heaviside(t-4999...
x(t) = 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*79991^(1/2)*Heaviside(t-4999...
x(t) = 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*79991^(1/2)*Heaviside(t-4999...
x(t) = 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*79991^(1/2)*Heaviside(t-4999...
x(t) = 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*79991^(1/2)*Heaviside(t-4999...
x(t) = 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*79991^(1/2)*Heaviside(t-4999...
x(t) = 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*79991^(1/2)*Heaviside(t-4999...
x(t) = 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*79991^(1/2)*Heaviside(t-4999...
 

このように超関数が入っていてもdsolveは解を与えます。 

右辺をunapplyを使って関数にします。 

> ans:=unapply(rhs(%),t);
 

proc (t) options operator, arrow; 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*7...
proc (t) options operator, arrow; 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*7...
proc (t) options operator, arrow; 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*7...
proc (t) options operator, arrow; 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*7...
proc (t) options operator, arrow; 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*7...
proc (t) options operator, arrow; 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*7...
proc (t) options operator, arrow; 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*7...
proc (t) options operator, arrow; 147/79991000*exp(-3/20*t)*sin(1/20*79991^(1/2)*t)*79991^(1/2)+49/1000*exp(-3/20*t)*cos(1/20*79991^(1/2)*t)-49/1000+1/10*Heaviside(t-4999999999/10000000000)-3/799910*7...
 

プロットしてみましょう。 

> plot({ans(t),g(t)},t=0..10);
 

Plot 

>
 

ステップ応答の様子がよくわかります。 

 

これをまたアニメーション表示してみましょう。 

> p:=[]:
 

pはアニメーションのフレームを追加していくための変数で、最初は空にしておきます。 

> r:=1:
 

rは車輪の半径です。 

> c:=cuboid([-0.2,-0.2,-0.2],[0.2,0.2,1.2],color=green):
 

cはタイヤにつながっている部分を直方体(cuboid)で作成しておきます。 

対角の2点の座標で定義します。色は緑(green)にします。 

これをtranslateでシフトしてアニメーションします。 

> M:=cuboid([-0.5,-0.5,0],[0.5,0.5,1],color=red):
 

Mは上に載っている物体を現す立方体です。色は赤にします。 

 

以下のforループで時間(tが0から1秒まで)ごとにアニメーションのフレームを作成します。 

スプリングはtubeplotを使って表現します。スプリングの長さをlで表します。 

hは車輪の中心の高さでHは上の質点(赤)の中心の高さです。 

> for t from 0 by 0.05 to 5 do
 h:=g(t):
 pp:=plot3d([x,y,g(x+t)],x=-Pi/2..Pi/2,y=-1..1,grid=[15,12],color=gray):
 pp:=pp,translate(rotate(cylinder([0,0,0],r,0.2,color=blue),Pi/2,t*2*Pi,0),0,-0.1,h+r):
 pp:=pp,translate(c,0,0,h+r):
 H:=ans(t):
 l:=H-h+1:
 tt:=tubeplot([cos(u*8*Pi)/4,sin(u*8*Pi)/4,u*l],u=0..1,radius=0.05,color=yellow):
 pp:=pp,translate(tt,0,0,h+2*r+0.2):
 pp:=pp,translate(M,0,0,h+2*r+0.2+l):
 p:=p,display(pp):
end:
 

できたフレームをアニメーション表示します。 

> display(p,insequence=true,scaling=constrained);
 

Plot