066_Vibration5.mw

機械振動論 

車体振動 

より具体的なシミュレーションにチャレンジしましょう。車で道路を走行するときの振動についてシミュレーションしてみましょう。 

アニメーション 

例によってアニメーションを作成してみましょう。 

 

まずMapleを初期化します。 

> restart:
 

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

> with(plots):with(plottools):
 

Warning, the name changecoords has been redefined 

Warning, the assigned name arrow now has a global binding 

道路の形状を表す関すを定義します。この場合正弦関数です。 

> g:=x->sin(x)/3;
 

proc (x) options operator, arrow; 1/3*sin(x) end proc 

車の速度、長さ、幅に関するパラメータを定義します。 

> V:=2*Pi;L:=2;W:=1.5;
 

2*Pi 

2 

1.5 

タイヤの半径のパラメータを決めています。 

> r:=1;
 

1 

タイヤを表すシリンダーを作っておきます。これを回転(rotate)と移動(translate)でアニメーションします。 

> ti1:=rotate(cylinder([0,0,0],r,0.4,color=blue),Pi/2,0,0):
 

車体を表す直方体です。これも回転と移動でアニメーションします。 

> body:=cuboid([-L-1,-W+0.5,-0.5],[L+1,W-0.5,0.5],color=green):
 

フレームを保存する変数pを空にセットしておきます。 

> p:=[];
 

[] 

時間(t)を進めながらフレームを作成していきます。 

道路をplot3dで表示します。 

h1は前輪の高さです。h2は後輪の高さです。 

車体の傾きはこのh1とh2から正接の逆関数(arctan)を使って計算しています。 

> for t from 0 by 0.02 to 1 do
 pp:=plot3d([u,v,g(V*t+u)],u=-6..6,v=-3..3):
 h1:=g(V*t+L):
 pp:=pp,translate(rotate(ti1,0,-2*t,0),L,W,h1+r):
 pp:=pp,translate(rotate(ti1,0,-2*t,0),L,-W,h1+r):
 h2:=g(V*t-L):
 pp:=pp,translate(rotate(ti1,0,-2*t,0),-L,W,h2+r):
 pp:=pp,translate(rotate(ti1,0,-2*t,0),-L,-W,h2+r):
 ang:=arctan(h1-h2,2*L):
 pp:=pp,translate(rotate(body,0,ang,0),0,0,2):
 p:=p,display(pp):
end:
 

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

Plot 

>
 

シミュレーション 

それでは実際にモデル化していきましょう。 

 

まず重力加速度のパラメータGを設定しておきます。 

> G:=9.8;
 

9.8 

変数tおよびVをリセットしておきます。 

> t:='t';V:='V';
 

t 

V 

車体の回転の運動方程式、前輪の運動方程式および後輪の運動方程式の三つの連立方程式でモデル化します。 

 

まずは車体の回転(ヨー)の運動方程式です。 

> eq1:=Iota*diff(phi(t),t$2)+c*diff(phi(t),t)+k*phi(t)=m1*diff(x1(t),t$2)-m2*diff(x2(t),t$2);
 

Iota*(diff(diff(phi(t), t), t))+c*(diff(phi(t), t))+k*phi(t) = m1*(diff(diff(x1(t), t), t))-m2*(diff(diff(x2(t), t), t))
Iota*(diff(diff(phi(t), t), t))+c*(diff(phi(t), t))+k*phi(t) = m1*(diff(diff(x1(t), t), t))-m2*(diff(diff(x2(t), t), t))
 

Iは車体の慣性モーメントです。phi(t)が時間tでの車体の前後の傾きの角度を単位ラジアンで表現します。 

これを一階微分(diff)したものが角速度です。更に微分(二階微分)したものが角加速度です。 

m1は車体の前側の質点の質量、m2は後ろ側の質量です。重心は前後の2質点の中心にあるものとしそこで回転の運動方程式が立てられます。 

前後の質点の力の差が車体を回転させようとする力となります。 

 

次が前輪の運動方程式です。 

x1(t)は前輪の高さです。車体の回転の力が加わります。 

> eq2:=m1*(diff(x1(t),t$2)+G)+c1*diff(x1(t),t)+k1*(x1(t)-g(V*t+L)+L*tan(phi(t)))=0;
 

m1*((diff(diff(x1(t), t), t))+9.8)+c1*(diff(x1(t), t))+k1*(x1(t)-1/3*sin(V*t+2)+2*tan(phi(t))) = 0
m1*((diff(diff(x1(t), t), t))+9.8)+c1*(diff(x1(t), t))+k1*(x1(t)-1/3*sin(V*t+2)+2*tan(phi(t))) = 0
 

同様に後輪の運動方程式です。 

> eq3:=m2*(diff(x2(t),t$2)+G)+c2*diff(x2(t),t)+k2*(x2(t)-g(V*t-L)-L*tan(phi(t)))=0;
 

m2*((diff(diff(x2(t), t), t))+9.8)+c2*(diff(x2(t), t))+k2*(x2(t)-1/3*sin(V*t-2)-2*tan(phi(t))) = 0
m2*((diff(diff(x2(t), t), t))+9.8)+c2*(diff(x2(t), t))+k2*(x2(t)-1/3*sin(V*t-2)-2*tan(phi(t))) = 0
 

速度Vを決めておきます。 

> V:=Pi;
 

Pi 

回転の慣性モーメント、粘性減衰定数、バネ定数を決めます。 

> Iota:=1;c:=0.001;k:=1;
 

1 

0.1e-2 

1 

車体の前側の質点の質量、粘性減衰定数、バネ定数の値を決めます。 

> m1:=1;c1:=0.1;k1:=100;
 

1 

.1 

100 

同じく後ろ側の質点の質量、粘性減衰定数、バネ定数の値を決めます。 

> m2:=1;c2:=0.1;k2:=100;
 

1 

.1 

100 

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

> eq1:=evalf(eq1);eq2:=evalf(eq2);eq3:=evalf(eq3);
 

(diff(diff(phi(t), t), t))+0.1e-2*(diff(phi(t), t))+phi(t) = (diff(diff(x1(t), t), t))-1.*(diff(diff(x2(t), t), t))
(diff(diff(phi(t), t), t))+0.1e-2*(diff(phi(t), t))+phi(t) = (diff(diff(x1(t), t), t))-1.*(diff(diff(x2(t), t), t))
 

(diff(diff(x1(t), t), t))+9.8+.1*(diff(x1(t), t))+100.*x1(t)-33.33333333*sin(3.141592654*t+2.)+200.*tan(phi(t)) = 0.
(diff(diff(x1(t), t), t))+9.8+.1*(diff(x1(t), t))+100.*x1(t)-33.33333333*sin(3.141592654*t+2.)+200.*tan(phi(t)) = 0.
(diff(diff(x1(t), t), t))+9.8+.1*(diff(x1(t), t))+100.*x1(t)-33.33333333*sin(3.141592654*t+2.)+200.*tan(phi(t)) = 0.
 

(diff(diff(x2(t), t), t))+9.8+.1*(diff(x2(t), t))+100.*x2(t)-33.33333333*sin(3.141592654*t-2.)-200.*tan(phi(t)) = 0.
(diff(diff(x2(t), t), t))+9.8+.1*(diff(x2(t), t))+100.*x2(t)-33.33333333*sin(3.141592654*t-2.)-200.*tan(phi(t)) = 0.
(diff(diff(x2(t), t), t))+9.8+.1*(diff(x2(t), t))+100.*x2(t)-33.33333333*sin(3.141592654*t-2.)-200.*tan(phi(t)) = 0.
 

初期値(角度、前輪の高さ、後輪の高さ)をきめます。 

> init:={phi(0)=0,x1(0)=0,x2(0)=0};
 

{x2(0) = 0, x1(0) = 0, phi(0) = 0} 

角速度の初期値、前輪の上下の速度および後輪の上下の速度の初期値をきめます。 

> init1:={D(phi)(0)=0,D(x1)(0)=0,D(x2)(0)=0};
 

{(D(x2))(0) = 0, (D(x1))(0) = 0, (D(phi))(0) = 0} 

計算のサンプル個数mと間隔dを定義します。 

> m:=50;d:=6*Pi/m;
 

50 

3/25*Pi 

サンプル時間のセットを作成します。 

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

準備ができましたのでdsolveで解きます。output=で計算されるサンプル時間sampを指定します。 

> ans:=dsolve({eq1,eq2,eq3} union init union init1,{phi(t),x1(t),x2(t)},numeric,output=samp):
 

結果をみていきましょう。 

計算されるデータのリストです。 

> seq(ans[1,1][j],j=1..7);
 

t, phi(t), diff(phi(t), t), x1(t), diff(x1(t), t), x2(t), diff(x2(t), t) 

角度の変化をプロットします。 

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

Plot 

前輪の高さの変化です。 

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

Plot 

同様に後輪の高さの変化です。 

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

Plot 

この結果をアニメーションしましょう。 

 

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

> p:=[];
 

[] 

時間を進めながらフレームを描きます。車体の角度はphi(t)を使用します。 

> for j from 0 to m do
 t:=ans[2,1][j+1,1]:
 pp:=plot3d([u,v,g(V*t+u)],u=-6..6,v=-3..3):
 h1:=g(V*t+L):
 pp:=pp,translate(rotate(ti1,0,-2*t,0),L,W,h1+r):
 pp:=pp,translate(rotate(ti1,0,-2*t,0),L,-W,h1+r):
 h2:=g(V*t-L):
 pp:=pp,translate(rotate(ti1,0,-2*t,0),-L,W,h2+r):
 pp:=pp,translate(rotate(ti1,0,-2*t,0),-L,-W,h2+r):
 h:=(ans[2,1][j+1,4]+ans[2,1][j+1,6])/2:
 ang:=ans[2,1][j+1,2]:
 pp:=pp,translate(rotate(body,0,ang,0),0,0,2+h):
 p:=p,display(pp):
end:
 

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

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

Plot 

>
 

もう一つのケース 

地面の形状を変えてシミュレーションを行ってみましょう。 

 

まずパラメータをリセットしておきます。 

> t:='t';g:='g';k:='k';k1:='k1';k2:='k2';L:='L';m1:='m1';m2:='m2';G:='G';c1:='c1';c2:='c2';V:='V';
 

t 

g 

k 

k1 

k2 

L 

m1 

m2 

G 

c1 

c2 

V 

車体回転の運動方程式です。 

> eq1:=Iota*diff(phi(t),t$2)+c*diff(phi(t),t)+k*phi(t)=m1*diff(x1(t),t$2)-m2*diff(x2(t),t$2);
 

(diff(diff(phi(t), t), t))+0.1e-2*(diff(phi(t), t))+k*phi(t) = m1*(diff(diff(x1(t), t), t))-m2*(diff(diff(x2(t), t), t))
(diff(diff(phi(t), t), t))+0.1e-2*(diff(phi(t), t))+k*phi(t) = m1*(diff(diff(x1(t), t), t))-m2*(diff(diff(x2(t), t), t))
 

前輪側の質点の運動方程式です。 

> eq2:=m1*(diff(x1(t),t$2)+G)+c1*diff(x1(t),t)+k1*(x1(t)-g(V*t+L)+L*tan(phi(t)))=0;
 

m1*((diff(diff(x1(t), t), t))+G)+c1*(diff(x1(t), t))+k1*(x1(t)-g(V*t+L)+L*tan(phi(t))) = 0
m1*((diff(diff(x1(t), t), t))+G)+c1*(diff(x1(t), t))+k1*(x1(t)-g(V*t+L)+L*tan(phi(t))) = 0
 

後輪側の質点の運動方程式です。 

> eq3:=m2*(diff(x2(t),t$2)+G)+c2*diff(x2(t),t)+k2*(x2(t)-g(V*t-L)-L*tan(phi(t)))=0;
 

m2*((diff(diff(x2(t), t), t))+G)+c2*(diff(x2(t), t))+k2*(x2(t)-g(V*t-L)-L*tan(phi(t))) = 0
m2*((diff(diff(x2(t), t), t))+G)+c2*(diff(x2(t), t))+k2*(x2(t)-g(V*t-L)-L*tan(phi(t))) = 0
 

速度のパラメータを設定します。 

> V:=Pi;
 

Pi 

重力加速度と車体の長さ(重心から質店までの距離)の設定をします。 

> G:=9.8;L:=2;
 

9.8 

2 

慣性モーメント、車体の回転の粘性減衰定数、バネ定数をセットします。 

> Iota:=1;c:=0.001;k:=100;
 

1 

0.1e-2 

100 

前側の質点の質量、粘性減衰定数およびバネ定数です。 

> m1:=1;c1:=0.1;k1:=1000;
 

1 

.1 

1000 

後ろ側の質点の質量、粘性減衰定数およびバネ定数です。 

> m2:=1;c2:=0.1;k2:=1000;
 

1 

.1 

1000 

サンプル個数と刻み(サンプル間隔)をセットします。 

> m:=100;d:=4*Pi/m;
 

100 

1/25*Pi 

サンプル時間のセットを生成します。 

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

道路のディップの場所のパラメータです。 

> LL:=4;
 

4 

道路の形状の関数です。区分関数(piecewise)を使っています。 

> g:=x->piecewise(x<LL,0,x<LL+0.5,(x-LL),x<LL+1,(-x+LL+1),0);
 

proc (x) options operator, arrow; piecewise(x < LL, 0, x < LL+.5, x-LL, x < LL+1, -x+LL+1, 0) end proc
proc (x) options operator, arrow; piecewise(x < LL, 0, x < LL+.5, x-LL, x < LL+1, -x+LL+1, 0) end proc
 

道路の形状をプロットします。 

> plot(g(x),x=0..m*d*V);
 

Plot 

数値化し方程式を簡単にします。 

> eq1:=evalf(eq1);eq2:=evalf(eq2);eq3:=evalf(eq3);
 

(diff(diff(phi(t), t), t))+0.1e-2*(diff(phi(t), t))+100.*phi(t) = (diff(diff(x1(t), t), t))-1.*(diff(diff(x2(t), t), t))
(diff(diff(phi(t), t), t))+0.1e-2*(diff(phi(t), t))+100.*phi(t) = (diff(diff(x1(t), t), t))-1.*(diff(diff(x2(t), t), t))
 

(diff(diff(x1(t), t), t))+9.8+.1*(diff(x1(t), t))+1000.*x1(t)-1000.*piecewise(3.141592654*t < 2., 0., 3.141592654*t < 2.5, 3.141592654*t-2., 3.141592654*t < 3., 3.-3.141592654*t, 0.)+2000.*tan(phi(t))...
(diff(diff(x1(t), t), t))+9.8+.1*(diff(x1(t), t))+1000.*x1(t)-1000.*piecewise(3.141592654*t < 2., 0., 3.141592654*t < 2.5, 3.141592654*t-2., 3.141592654*t < 3., 3.-3.141592654*t, 0.)+2000.*tan(phi(t))...
(diff(diff(x1(t), t), t))+9.8+.1*(diff(x1(t), t))+1000.*x1(t)-1000.*piecewise(3.141592654*t < 2., 0., 3.141592654*t < 2.5, 3.141592654*t-2., 3.141592654*t < 3., 3.-3.141592654*t, 0.)+2000.*tan(phi(t))...
(diff(diff(x1(t), t), t))+9.8+.1*(diff(x1(t), t))+1000.*x1(t)-1000.*piecewise(3.141592654*t < 2., 0., 3.141592654*t < 2.5, 3.141592654*t-2., 3.141592654*t < 3., 3.-3.141592654*t, 0.)+2000.*tan(phi(t))...
 

(diff(diff(x2(t), t), t))+9.8+.1*(diff(x2(t), t))+1000.*x2(t)-1000.*piecewise(3.141592654*t < 6., 0., 3.141592654*t < 6.5, 3.141592654*t-6., 3.141592654*t < 7., 7.-3.141592654*t, 0.)-2000.*tan(phi(t))...
(diff(diff(x2(t), t), t))+9.8+.1*(diff(x2(t), t))+1000.*x2(t)-1000.*piecewise(3.141592654*t < 6., 0., 3.141592654*t < 6.5, 3.141592654*t-6., 3.141592654*t < 7., 7.-3.141592654*t, 0.)-2000.*tan(phi(t))...
(diff(diff(x2(t), t), t))+9.8+.1*(diff(x2(t), t))+1000.*x2(t)-1000.*piecewise(3.141592654*t < 6., 0., 3.141592654*t < 6.5, 3.141592654*t-6., 3.141592654*t < 7., 7.-3.141592654*t, 0.)-2000.*tan(phi(t))...
(diff(diff(x2(t), t), t))+9.8+.1*(diff(x2(t), t))+1000.*x2(t)-1000.*piecewise(3.141592654*t < 6., 0., 3.141592654*t < 6.5, 3.141592654*t-6., 3.141592654*t < 7., 7.-3.141592654*t, 0.)-2000.*tan(phi(t))...
 

dsolveで方程式を解きます。 

> ans:=dsolve({eq1,eq2,eq3} union init union init1,{phi(t),x1(t),x2(t)},numeric,output=samp):
 

車体の傾きの角度の変化です。 

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

Plot 

前輪側の質点の高さの変化です。 

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

Plot 

後輪側の質店の高さのプロットです。 

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

Plot 

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

> p:=[];
 

[] 

時間を進めながらフレームを作成します。 

> for j from 0 to m do
 t:=ans[2,1][j+1,1]:
 pp:=plot3d([u,v,g(V*t+u)],u=-6..6,v=-3..3):
 h1:=g(V*t+L):
 pp:=pp,translate(rotate(ti1,0,-2*t,0),L,W,h1+r):
 pp:=pp,translate(rotate(ti1,0,-2*t,0),L,-W,h1+r):
 h2:=g(V*t-L):
 pp:=pp,translate(rotate(ti1,0,-2*t,0),-L,W,h2+r):
 pp:=pp,translate(rotate(ti1,0,-2*t,0),-L,-W,h2+r):
 h:=(ans[2,1][j+1,4]+ans[2,1][j+1,6])/2:
 ang:=ans[2,1][j+1,2]:
 pp:=pp,translate(rotate(body,0,ang,0),0,0,2+h):
 p:=p,display(pp):
end:
 

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

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

Plot