068_Vibration7.mw

機械振動論 

固有振動と振動モード 

物体の振動は振動モードと呼ばれる特定のパターンで振動します。 

実際の振動は複数の振動モードが組み合わされています。 

板バネの振動 

まずは板バネの振動のシミュレーションを行ってみましょう。 

 

Mapleを初期化します。 

> restart:
 

板バネを30分割してモデル化します。 

> n:=30;
 

30

分割したそれぞれの板に関して運動方程式を生成します。 

一方の端は固定しないため最後の板を除き同じパターンの運動方程式です。 

 

x[j]は各板の高さを表します。一階微分は速度、二階微分は加速度です。 

> for j to n-1 do
 eq[j]:=  m*diff(x[j](t),t$2)+k*(x[j](t)-x[j-1](t))-k*(x[j+1](t)-x[j](t));
end:
 

端(30番目)の板は隣がありませんので別に定義します。 

> eq[30]:=m*diff(x[30](t),t$2)+k*(x[30](t)-x[29](t));
 

m*(diff(diff(x[30](t), t), t))+k*(x[30](t)-x[29](t)) 

各板の質量mを決めます。 

> m:=500;
 

500 

同じくバネ定数を決めます。 

> k:=1.0e7;
 

0.10e8 

板バネの元(x[0])は固定にします。x[0]は常に0です。 

> x[0]:=t->0;
 

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

初期値の設定をします。まず空集合にしておきます。 

> init:={};
 

{} 

30番目(端)以外はすべて初期高さ(x[j])および初速(D(x[j])は全て0にします。 

> for j to n-1 do
 init:=init union {x[j](0)=0,D(x[j])(0)=0}
end:
 

30番目の板は初期高さは0で初速を0.1(m/s:単位系MKS)とします。 

> init:=init union {x[30](0)=0,D(x[30])(0)=0.1}:
 

方程式の集合をセットします。 

> equ:={seq(eq[j],j=1..n)}:
 

サンプリングの個数mmを決めます。 

> mm:=50;
 

50 

シミュレーション時間TTを決めます。 

> TT:=1;
 

1 

出力刻みが決まります。 

> d:=TT/mm;
 

1/50 

サンプリングのセットを生成します。 

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

dsolveで運動方程式の連立微分方程式系を解きます。 

> ans:=dsolve(equ union init,{seq(x[j](t),j=1..n)},numeric,output=samp):
 

求められる解は次のものです。 

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

t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t), diff(x[7](t), t), x[8](t),...
t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t), diff(x[7](t), t), x[8](t),...
t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t), diff(x[7](t), t), x[8](t),...
t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t), diff(x[7](t), t), x[8](t),...
t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t), diff(x[7](t), t), x[8](t),...
t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t), diff(x[7](t), t), x[8](t),...
t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t), diff(x[7](t), t), x[8](t),...
t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t), diff(x[7](t), t), x[8](t),...
t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t), diff(x[7](t), t), x[8](t),...
t, x[1](t), diff(x[1](t), t), x[2](t), diff(x[2](t), t), x[3](t), diff(x[3](t), t), x[4](t), diff(x[4](t), t), x[5](t), diff(x[5](t), t), x[6](t), diff(x[6](t), t), x[7](t), diff(x[7](t), t), x[8](t),...
 

結果をプロットしましょう。 

> with(plots):with(plottools):
 

Warning, the name changecoords has been redefined 

Warning, the assigned name arrow now has a global binding 

解(ans)の二番目はx[1]の時間に対する変化を表します。 

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

Plot 

60番目はx[30]すなわち一番端(フリー)の時間に対する変化です。 

> listplot([seq(ans[2,1][j,60],j=1..mm+1)]);
 

Plot 

板全体の動きをアニメーションしましょう。 

 

各板のオブジェクトを直方体で用意します。 

> cb:=cuboid([0,0,0],[4,1,0.1]):
 

表示してみます。 

> display(cb,scaling=constrained);
 

Plot 

> kabe:=cuboid([0,-1,-3],[4,0,3],color=gray):
 

フレームを保存する変数pを初期化します。 

> p:=[];
 

[] 

時間変化の各フレームを作成します。変位は強調(1500倍)してあります。 

> for j to mm+1 do
 pp:=[kabe,cb]:
 for i to n do
   pp:=pp,translate(cb,0,i,1500*ans[2,1][j,2*i]):
 end:
 p:=p,display(pp):
end:
 

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

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

Plot 

>
 

固有値問題 

この板バネの系の固有値問題を解き固有モードと振動スペクトラムを求めます。 

 

これらの30個の運動方程式は行列形式で表現することができます。 

 

[M]*{diff(x, x, x)}+[K]*{x} = {f} 

 

Mを質量行列、Kを剛性行列と呼びます。 xを変位ベクトル、fを力ベクトルと呼びます。 

 

このときすべての自由度が同位相で振動するとすると 

 

Omega^2*[M]*X-[K]*X = 0 ([K]*X = lambda*[M]*X) 

 

という行列方程式となります。これを一般固有値問題(A*x = lambda*B*x)とよんでいます。ここでOmegaは固有角振動数です。 

 

[M]の逆行列を両辺にかければ通常の標準固有値問題(A*x = lambda*x)となります。 

 

固有値問題を解く為に線形代数パッケージを呼び出します。 

 

> with(LinearAlgebra):
 

剛体行列を生成します。対角が質量の単純な対角行列です。 

> M:=DiagonalMatrix([seq(m,j=1..n)]);
 

Matrix(%id = 680621340) 

剛性行列を生成します。まず対角要素を入れます。 

> K:=DiagonalMatrix([seq(2*k,j=1..n)]);
 

Matrix(%id = 681069184) 

上側と下側のsubdiagonal(対角のひとつ上、ひとつ下)を設定します。 

> for j to n-1 do
 K[j,j+1]:=-k:
 K[j+1,j]:=-k:
end:
 

行列KとMを与えて固有値(W)と固有ベクトル(P)を関数Eigenvectorsで求めます。 

> (P,W):=Eigenvectors(K,M,output=['vectors', 'values']);
 

Matrix(%id = 681982740), Vector[column](%id = 679119820) 

求められた固有値です。 

> seq([j,W[j]],j=1..n);
 

[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
[1, 818.802349900219952+0.*I], [2, 205.227064324200000+0.*I], [3, 1834.42974399800005+0.*I], [4, 3241.68753519080019+0.*I], [5, 5026.13535421679990+0.*I], [6, 7169.46235170900036+0.*I], [7, 9649.67509...
 

複素数になっていますが虚部はすべて0です。 

 

固有ベクトル(実部)をすべてプロットしてみましょう。 

この行列は固有モード行列あるいはモード行列と呼ばれます。 

> matrixplot(map(Re,P),shading=zhue);
 

Plot 

周波数がだんだん高くなっています。 

 

最初の振動モード(最も低周波)をプロットしてみましょう。 

このモードは一次モードあるいは基本モードとも呼ばれます。 

> listplot([seq(Re(P[j,1]),j=1..n)]);
 

Plot 

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

まずはフレームを保存する変数を空にしておきます。 

> p:=[];
 

[] 

時間変化の各フレームを生成します。 

> for j from 0 to mm  do
 pp:=[kabe,cb]:
 for i to n do
   pp:=pp,translate(cb,0,i,10*Re(P[i,1])*sin(2*Pi*j/mm)):
 end:
 p:=p,display(pp):
end:
 

アニメーション表示します。 

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

Plot 

2番目の振動モードをプロットします。 

> listplot([seq(Re(P[j,4]),j=1..n)]);
 

Plot 

アニメーション表示します。 

> p:=[];
 

[] 

> for j from 0 to mm  do
 pp:=[kabe,cb]:
 for i to n do
   pp:=pp,translate(cb,0,i,10*Re(P[i,4])*sin(2*Pi*j/mm)):
 end:
 p:=p,display(pp):
end:
 

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

Plot 

3番目の振動モードです。 

> listplot([seq(Re(P[j,6]),j=1..n)]);
 

Plot 

> p:=[];
 

[] 

> for j from 0 to mm  do
 pp:=[kabe,cb]:
 for i to n do
   pp:=pp,translate(cb,0,i,6*Re(P[i,6])*sin(2*Pi*j/mm)):
 end:
 p:=p,display(pp):
end:
 

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

Plot 

4番目の振動モードです。 

> listplot([seq(Re(P[j,8]),j=1..n)]);
 

Plot 

> p:=[];
 

[] 

> for j from 0 to mm  do
 pp:=[kabe,cb]:
 for i to n do
   pp:=pp,translate(cb,0,i,3*Re(P[i,8])*sin(2*Pi*j/mm)):
 end:
 p:=p,display(pp):
end:
 

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

Plot 

5番目の振動モードです。 

> listplot([seq(Re(P[j,10]),j=1..n)]);
 

Plot 

> p:=[];
 

[] 

> for j from 0 to mm  do
 pp:=[kabe,cb]:
 for i to n do
   pp:=pp,translate(cb,0,i,3*Re(P[i,10])*sin(2*Pi*j/mm)):
 end:
 p:=p,display(pp):
end:
 

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

Plot 

2番目と3番目を一緒にアニメーションしてみます。 

> p:=[];
 

[] 

> for j from 0 to mm  do
 pp:=[kabe,cb]:
 pp:=pp,translate(kabe,5,0,0):
 pp:=pp,translate(cb,5,0,0):
 for i to n do
   pp:=pp,translate(cb,0,i,6*Re(P[i,4])*sin(2*Pi*j/mm)):
   pp:=pp,translate(cb,5,i,6*Re(P[i,6])*sin(2*Pi*j/mm)):
 end:
 p:=p,display(pp):
end:
 

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

Plot 

各モードの固有角振動(Omega)をプロットしてみます(スペクトラム)。 

> listplot([seq(sqrt(Re(W[j])),j=1..n)]);
 

Plot