065_Vibration4.mw

機械振動論 

ねじり運動 

振動には物体のねじれに対応する振動があります。 

ここではねじれ振動について考えてみましょう。 

アニメーション 

視覚的に理解するためにアニメーションを作成します。 

 

まずMapleを初期化します。 

> restart:
 

プロット用のパッケージを呼び出しておきます。 

> with(plots):with(plottools):
 

Warning, the name changecoords has been redefined 

Warning, the assigned name arrow now has a global binding 

スプリングをつっている天井のオブジェクトを作成します。 

> c:=cuboid([-3,-1,4],[-5,1,5]):
 

強制振動(回転)する物体を直方体で用意します。 

> c0:=cuboid([-5,-6,-1],[5,-2,1]):
 

それにつながれた解析対象の回転物体のオブジェクトを用意します。 

> c1:=cuboid([-5,-1,-1],[5,1,1]):
 

スプリングの長さをきめます。 

> ll:=3;
 

3 

スプリングを作成します。 

> sp:=tubeplot([cos(5*Pi*u)-4,sin(5*Pi*u),u*ll+1],u=0..1,radius=0.2):
 

回転の中心の棒を用意します。 

> bou:=translate(rotate(cylinder([0,0,0],0.2,6),Pi/2,0,0),0,-3,0):
 

針の先の変化の面を作成します。 

> hei:=plot3d([u,v,sin(v)/2-ll-1],u=3..6,v=-2*Pi..2*Pi,grid=[10,20]):
 

針をコーンで定義します。 

> kone:=cone([4,0,-4],1,ll):
 

表示してみましょう。 

> display({c0,c1,c,bou,sp,hei,kone},scaling=constrained);
 

Plot 

入力変化面を動かし物体が回転するアニメーションを作成します。 

 

アニメーションのフレームを保存していく変数を空で用意します。 

> p:=[];
 

[] 

時間tを動かしながらフレームを作成していきます。 

回転はrotateを使って行います。 

> for t from 0 by 0.2 to evalf(2*Pi) do
 pp:=c0:
 pp:=pp,bou:
 pp:=pp,plot3d([u,v,sin(v+t)/2-ll-1],u=3..6,v=-2*Pi..2*Pi,grid=[10,20]):
 h:=evalf(sin(t)/2):
 pp:=pp,translate(kone,0,0,h):
 phi:=arctan(h,4):
 pp:=pp,rotate(c1,0,phi,0):
 pp:=pp,rotate(c,0,1.5*phi,0):
 l:=ll+h:
 sp:=tubeplot([cos(5*Pi*u)-4,sin(5*Pi*u),u*l+1],u=0..1,radius=0.2):
 pp:=pp,translate(sp,0,0,-h):
 p:=p,display(pp):
end:
 

変化面はh:=sin(t)/2で変化します。 

φは回転角(単位:ラジアン)です。正接(tan)の逆関数で計算しています。 

 

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

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

Plot 

>
 

ねじれ(回転)の運動方程式(自由度1) 

Mapleを初期化します。 

> restart:
 

回転の運動方程式もいままでの運動方程式と類似的な式となります。 

> Iota*diff(phi(t),t$2)+c*diff(phi(t),t)+k*phi(t)=0;
 

Iota*(diff(diff(phi(t), t), t))+c*(diff(phi(t), t))+k*phi(t) = 0 

位置の変わりに回転角に関する微分方程式となります。 

lは慣性モーメントで通常の運動方程式の質量Mに対応します。 

角度の二階微分は角加速度と呼ばれます。 

一階微分は角速度です。すなわち回転の速さです。 

粘性減衰定数cは回転に対する減衰を表します。 

kも同様に回転に関するバネ定数です。回転を復元しようとする力を表現します。 

 

同様にMapleではdsolveでとくことができます。 

> dsolve(%,phi(t));
 

phi(t) = _C1*exp(-1/2*(c-(c^2-4*k*Iota)^(1/2))*t/Iota)+_C2*exp(-1/2*(c+(c^2-4*k*Iota)^(1/2))*t/Iota) 

ここで_C1および_C2は任意定数で解は二つの解の線形結合となっています。 

>
 

多段のねじり振動 

次にねじり振動の多段のモデルを考えてみましょう。 

 

まずMapleを初期化します。 

> restart:
 

プロット用のパッケージを呼び出しておきます。 

> with(plots):with(plottools):
 

Warning, the name changecoords has been redefined 

Warning, the assigned name arrow now has a global binding 

段数を指定する変数nをセットします。 

今回は三段です。 

> n:=3;
 

3 

回転の運動方程式を生成します。 

> for j from 1 to n do
 eq[j]:=Iota[j]*diff(phi[j](t),t$2)+c[j]*diff(phi[j](t),t)+k[j]*(phi[j](t)-phi[j-1](t))=0;
end;
 

Iota[1]*(diff(diff(phi[1](t), t), t))+c[1]*(diff(phi[1](t), t))+k[1]*(phi[1](t)-phi[0](t)) = 0 

Iota[2]*(diff(diff(phi[2](t), t), t))+c[2]*(diff(phi[2](t), t))+k[2]*(phi[2](t)-phi[1](t)) = 0 

Iota[3]*(diff(diff(phi[3](t), t), t))+c[3]*(diff(phi[3](t), t))+k[3]*(phi[3](t)-phi[2](t)) = 0 

各定数をきめておきます。まず慣性モーメントです。 

> for j to n do
 Iota[j]:=1:
end:
 

次は粘性減衰定数です。 

> for j to n do
 c[j]:=0.3:
end:
 

同じくバネ定数です。 

> for j to n do
 k[j]:=2000:
end:
 

入力振動(強制振動)の関数です。 

> phi[0]:=t->sin(t)/5;
 

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

evalfで数値化して簡単化します。 

> for j to n do
 eq[j]:=evalf(eq[j]);
end;
 

(diff(diff(phi[1](t), t), t))+.3*(diff(phi[1](t), t))+2000.*phi[1](t)-400.*sin(t) = 0.
(diff(diff(phi[1](t), t), t))+.3*(diff(phi[1](t), t))+2000.*phi[1](t)-400.*sin(t) = 0.
 

(diff(diff(phi[2](t), t), t))+.3*(diff(phi[2](t), t))+2000.*phi[2](t)-2000.*phi[1](t) = 0.
(diff(diff(phi[2](t), t), t))+.3*(diff(phi[2](t), t))+2000.*phi[2](t)-2000.*phi[1](t) = 0.
 

(diff(diff(phi[3](t), t), t))+.3*(diff(phi[3](t), t))+2000.*phi[3](t)-2000.*phi[2](t) = 0.
(diff(diff(phi[3](t), t), t))+.3*(diff(phi[3](t), t))+2000.*phi[3](t)-2000.*phi[2](t) = 0.
 

初期値を設定します。 

各オブジェクトの回転角の初期値です。すべて0にします。 

> init:={}:
 

> for j to n do
 init:=init union {phi[j](0)=0}
end:
 

初期の角速度です。同じく0にします。 

> init1:={}:
 

> for j to n do
 init1:=init1 union {D(phi[j])(0)=0}
end:
 

連立方程式の集合です。 

> equ:={}:
 

> for j to n do
 equ:=equ union {eq[j]}
end:
 

解くべき関数(回転角度)の集合です。 

> funcs:={}:
 

> for j to n do
 

> funcs:=funcs union {phi[j](t)}
 

> end:
 

計算個数と計算ステップ幅をきめます。 

> m:=50;d:=evalf(2*Pi/m);
 

50 

.1256637062 

入力関数をプロットしておきましょう。 

> plot(phi[0](t),t=0..d*m);
 

Plot 

計算サンプル時間を生成します。dsolveのoutput=で指定します。 

> samp:=array([seq(d*j,j=0..m)]):
 

方程式をdsolveで解きます。numericを指定し数値解を得ます。outputで解くべきサンプル時間を指定します。 

> ans:=dsolve(equ union init union init1,funcs,numeric,output=samp):
 

出力されたデータのリストです。時間および角度と角速度です。 

> seq(ans[1,1][j],j=1..2*n+1);
 

t, phi[1](t), diff(phi[1](t), t), phi[2](t), diff(phi[2](t), t), phi[3](t), diff(phi[3](t), t) 

一段目の角度の時間変化のプロットをlistplotで行います。 

> listplot([seq(ans[2,1][j,2],j=1..m)]);
 

Plot 

同様に二段目の角度の時間変化です。 

> listplot([seq(ans[2,1][j,4],j=1..m)]);
 

Plot 

同じく3段目です。 

> listplot([seq(ans[2,1][j,6],j=1..m)]);
 

Plot 

>
 

このシミュレーション結果を元にアニメーションを作成しましょう。 

 

スプリングをつっている天井のオブジェクトを作成します。 

> c:=cuboid([-3,-1,4],[-5,1,5]):
 

回転する物体を用意します。 

> c0:=cuboid([-5,-1,-1],[5,1,1]):
 

> c1:=cuboid([-5,-2,-1],[5,2,1]):
 

> c2:=cuboid([-5,-2,-1],[5,2,1]):
 

> c3:=cuboid([-5,-2,-1],[5,2,1]):
 

スプリングの長さをきめます。 

> ll:=3;
 

3 

スプリングを作成します。 

> sp:=tubeplot([cos(5*Pi*u)-4,sin(5*Pi*u),u*ll+1],u=0..1,radius=0.2):
 

回転の中心の棒を用意します。 

> bou:=translate(rotate(cylinder([0,0,0],0.2,25),Pi/2,0,0),0,-20,0):
 

針の先の変化の面を作成します。 

> hei:=plot3d([u,v,sin(v)/2-ll-1],u=3..6,v=-2*Pi..2*Pi,grid=[10,20]):
 

針をコーンで定義します。 

> kone:=cone([4,0,-4],1,ll):
 

表示してみましょう。 

> display({c,c0,translate(c1,0,-6,0),translate(c2,0,-12,0),translate(c3,0,-18,0),bou,sp,hei,kone},scaling=constrained);
 

Plot 

入力変化面を動かし物体が回転するアニメーションを作成します。 

 

アニメーションのフレームを保存していく変数を空で用意します。 

> p:=[];
 

[] 

時間tを動かしながらフレームを作成していきます。 

> for j from 1 to m do
 pp:=c:
 pp:=pp,bou:
 t:=ans[2,1][j,1]:
 pp:=pp,plot3d([u,v,4*phi[0](t+v)-ll-1],u=3..6,v=-2*Pi..2*Pi,grid=[10,20]):
 h:=4*phi[0](t):
 pp:=pp,translate(kone,0,0,h):
 pp:=pp,rotate(c0,0,phi[0](t),0):
 phi1:=ans[2,1][j,2]:phi2:=ans[2,1][j,4]:phi3:=ans[2,1][j,6]:
 pp:=pp,translate(rotate(c1,0,phi1,0),0,-6,0):
 pp:=pp,translate(rotate(c2,0,phi2,0),0,-12,0):
 pp:=pp,translate(rotate(c3,0,phi3,0),0,-18,0):
 l:=ll+h:
 sp:=tubeplot([cos(5*Pi*u)-4,sin(5*Pi*u),u*l+1],u=0..1,radius=0.2):
 pp:=pp,translate(sp,0,0,-h):
 p:=p,display(pp):
end:
 

変化面はh:=sin(t)/2で変化します。 

φは回転角(単位:ラジアン)です。正接(tan)の逆関数で計算しています。 

 

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

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

Plot