034_DampingProblem.mw

減衰問題

イントロダクション

大きな回転する質量をもつ機械システムを考えます。 これが弾性物質の上に置かれているとします。減衰振幅を、周波数と減衰の関数として可視化することにより、 減衰 c の特性を決定することにします。

[Inserted Image]

> restart;

システムのパラメータ

機械の重さは m[M] = 1000kg、回転質量の重さは m[U] = 500kg です。弾性物質のばね定数は d = 100000 [ N/m ]。クランクシャフトの半径を r とし、その角速度を omega*t とします。

> val:= m[M]=1000, m[U]=500, d=100000:

機械の振動は、次の微分方程式で表されます。ここで、c は減衰係数を表します。

> deq := (m[M]+m[U])*diff(x(t),t,t) + c*diff(x(t),t) + d*x(t) + m[U]*r*omega^2*cos(omega*t) = 0:

この微分方程式の解を求めます。

> sol:=dsolve(deq,x(t));

sol := x(t) = exp(-1/2*(c-(c^2-4*d*m[M]-4*d*m[U])^(1/2))*t/(m[M]+m[U]))*_C2+exp(-1/2*(c+(c^2-4*d*m[M]-4*d*m[U])^(1/2))*t/(m[M]+m[U]))*_C1-((-omega^2*m[M]+d-omega^2*m[U])*cos(omega*t)+omega*sin(omega*t...sol := x(t) = exp(-1/2*(c-(c^2-4*d*m[M]-4*d*m[U])^(1/2))*t/(m[M]+m[U]))*_C2+exp(-1/2*(c+(c^2-4*d*m[M]-4*d*m[U])^(1/2))*t/(m[M]+m[U]))*_C1-((-omega^2*m[M]+d-omega^2*m[U])*cos(omega*t)+omega*sin(omega*t...

私たちは静止した状態だけに興味があるので、指数関数(過渡影響を表す)は無視することができます。

> x_stat:=unapply(subs(_C1=0,_C2=0,op(2,sol)),t);

x_stat := proc (t) options operator, arrow; -((-omega^2*m[M]+d-omega^2*m[U])*cos(omega*t)+omega*sin(omega*t)*c)*m[U]*r*omega^2/(omega^4*m[U]^2+(2*omega^4*m[M]-2*d*omega^2)*m[U]+omega^4*m[M]^2-2*d*m[M]...

土台に与える力は

> Fu := d*x_stat(t)+c*diff(x_stat(t),t);

Fu := -d*((-omega^2*m[M]+d-omega^2*m[U])*cos(omega*t)+omega*sin(omega*t)*c)*m[U]*r*omega^2/(omega^4*m[U]^2+(2*omega^4*m[M]-2*d*omega^2)*m[U]+omega^4*m[M]^2-2*d*m[M]*omega^2+d^2+c^2*omega^2)-c*(-(-omeg...Fu := -d*((-omega^2*m[M]+d-omega^2*m[U])*cos(omega*t)+omega*sin(omega*t)*c)*m[U]*r*omega^2/(omega^4*m[U]^2+(2*omega^4*m[M]-2*d*omega^2)*m[U]+omega^4*m[M]^2-2*d*m[M]*omega^2+d^2+c^2*omega^2)-c*(-(-omeg...

振幅を決めるため、次の加法定理を使い、振動する部分をその極大値に設定します。

> a*sin(omega*t)+b*cos(omega*t)=sqrt(a^2+b^2)*sin(omega*t+arcsin(b/sqrt(a^2+b^2)));

a*sin(omega*t)+b*cos(omega*t) = (a^2+b^2)^(1/2)*sin(omega*t+arcsin(b/(a^2+b^2)^(1/2)))

> match(subs(omega*t=ttt,Fu)=a*sin(ttt)+b*cos(ttt),ttt,'coef');

true

> coef;

{a = -c*omega^5*m[U]*r*(m[M]+m[U])/(omega^4*m[M]^2-2*d*m[M]*omega^2+2*m[M]*omega^4*m[U]+d^2-2*d*omega^2*m[U]+c^2*omega^2+omega^4*m[U]^2), b = omega^2*r*m[U]*(d*m[M]*omega^2-d^2+d*omega^2*m[U]-c^2*omeg...{a = -c*omega^5*m[U]*r*(m[M]+m[U])/(omega^4*m[M]^2-2*d*m[M]*omega^2+2*m[M]*omega^4*m[U]+d^2-2*d*omega^2*m[U]+c^2*omega^2+omega^4*m[U]^2), b = omega^2*r*m[U]*(d*m[M]*omega^2-d^2+d*omega^2*m[U]-c^2*omeg...

これにより、減衰振幅が与えられます。

> Af := unapply(simplify(subs(coef, sqrt(a^2+b^2)), symbolic), omega);

Af := proc (omega) options operator, arrow; omega^2*m[U]*r*(d^2+c^2*omega^2)^(1/2)/(omega^4*m[M]^2-2*d*m[M]*omega^2+2*m[M]*omega^4*m[U]+d^2-2*d*omega^2*m[U]+c^2*omega^2+omega^4*m[U]^2)^(1/2) end proc

(減衰なしの)アンバランスによる力は

> Fu := t -> m[U]*r*omega^2*cos(omega*t);

Fu := proc (t) options operator, arrow; m[U]*r*omega^2*cos(omega*t) end proc

よって、アンバランスによる振幅は

> Au := unapply(subs(cos(omega*t)=1, Fu(t)), omega);

Au := proc (omega) options operator, arrow; omega^2*m[U]*r end proc

減衰の特性は、振幅の比, Af/Au により、決まります。

> A := unapply(normal(Af(omega)/Au(omega)), omega, c);

A := proc (omega, c) options operator, arrow; (d^2+c^2*omega^2)^(1/2)/(omega^4*m[M]^2-2*d*m[M]*omega^2+2*m[M]*omega^4*m[U]+d^2-2*d*omega^2*m[U]+c^2*omega^2+omega^4*m[U]^2)^(1/2) end proc

c に対する、最大の振幅を求めます。

> tmp1 := [solve(diff(A(omega,c),omega)=0, omega)];

tmp1 := [0, (-(m[M]+m[U])*(d*m[M]+d*m[U]-(d^2*m[M]^2+2*d^2*m[M]*m[U]+d^2*m[U]^2+2*d*m[U]*c^2+2*d*m[M]*c^2)^(1/2))*d)^(1/2)/(c*(m[M]+m[U])), -(-(m[M]+m[U])*(d*m[M]+d*m[U]-(d^2*m[M]^2+2*d^2*m[M]*m[U]+d^...tmp1 := [0, (-(m[M]+m[U])*(d*m[M]+d*m[U]-(d^2*m[M]^2+2*d^2*m[M]*m[U]+d^2*m[U]^2+2*d*m[U]*c^2+2*d*m[M]*c^2)^(1/2))*d)^(1/2)/(c*(m[M]+m[U])), -(-(m[M]+m[U])*(d*m[M]+d*m[U]-(d^2*m[M]^2+2*d^2*m[M]*m[U]+d^...tmp1 := [0, (-(m[M]+m[U])*(d*m[M]+d*m[U]-(d^2*m[M]^2+2*d^2*m[M]*m[U]+d^2*m[U]^2+2*d*m[U]*c^2+2*d*m[M]*c^2)^(1/2))*d)^(1/2)/(c*(m[M]+m[U])), -(-(m[M]+m[U])*(d*m[M]+d*m[U]-(d^2*m[M]^2+2*d^2*m[M]*m[U]+d^...tmp1 := [0, (-(m[M]+m[U])*(d*m[M]+d*m[U]-(d^2*m[M]^2+2*d^2*m[M]*m[U]+d^2*m[U]^2+2*d*m[U]*c^2+2*d*m[M]*c^2)^(1/2))*d)^(1/2)/(c*(m[M]+m[U])), -(-(m[M]+m[U])*(d*m[M]+d*m[U]-(d^2*m[M]^2+2*d^2*m[M]*m[U]+d^...

物理的な現実性を考慮し、正の根を選びます。

> tmp2:=select(x->evalb(Re(subs(val,c=1.0,x))>0),tmp1)[1];

tmp2 := (-(m[M]+m[U])*(d*m[M]+d*m[U]-(d^2*m[M]^2+2*d^2*m[M]*m[U]+d^2*m[U]^2+2*d*m[U]*c^2+2*d*m[M]*c^2)^(1/2))*d)^(1/2)/(c*(m[M]+m[U]))

> tmp:=unapply(subs(val,tmp2),c);

tmp := proc (c) options operator, arrow; 1/1500*(-22500000000000000+150000000*(22500000000000000+300000000*c^2)^(1/2))^(1/2)/c end proc

> Max_A:=unapply(simplify(subs(val,A(tmp(c),c))),c);

Max_A := proc (c) options operator, arrow; (225000000+3*c^2)^(1/4)/((675000000000000000000-45000000000000000*(225000000+3*c^2)^(1/2)+9000000000000*c^2-300000000*c^2*(225000000+3*c^2)^(1/2)+(225000000+...Max_A := proc (c) options operator, arrow; (225000000+3*c^2)^(1/4)/((675000000000000000000-45000000000000000*(225000000+3*c^2)^(1/2)+9000000000000*c^2-300000000*c^2*(225000000+3*c^2)^(1/2)+(225000000+...

解のグラフィック表示

振幅の式にパラメータを代入し、それをグラフィック表示します。

> AA:=unapply(subs(val,A(omega,c)),omega,c):

> pl1 := plot3d(AA(omega,c), omega=0..20, c=2000..20000, colour=[(AA(omega,c)),0,-(AA(omega,c))]):
pl2 := plot3d([tmp(c),c,Max_A(c)], c=2000..20000, dummy=0..1, colour=green, thickness=3):
plots[display](pl1,pl2, labels=["omega [s]","c [kg/s]","Af/Au"], title="Relative Amplitude Vs. Frequency and Damping", axes=BOXED);

[Plot]

結論

Maple は、微分方程式の厳密解を与えてくれます。このシンボリックな解を使い、振動の振幅と土台に伝わる振幅の比を、計算したり、可視化したりできます。 この比は減衰の特性を知るにとても重要です。

参考文献

Heinz Waller and Reinhard Schmidt, Schwingungslehre f?r Ingenieure, Wissenschaftsverlag.