035_OscillationsBridge.mw

ローディングブリッジの振動

イントロダクション

このワークシートでは、ローディングブリッジに、質量 m を積んだ車両がのっている状態を考えます。質量 m が突然 離れたときの、ローディングブリッジの振動を決定します。

[Inserted Image]

> restart;

問題の定義

橋の長さは L、車両の外側の車輪の距離は l です。橋の左端から、車両中央までの距離は x0 で、その質量は mv です。E は弾性率、I はローディングブリッジの 断面2次モーメントです。
パラメータの単位は次の通りです。 長さ L, l, x0 は [ m ]、質量 m, mv は [ kg ] 、EI は [ Nm2 ] です。

パラメータの値を以下に示します。

> L:=10:
val:= l=1, x0=8, mv=1000, m=4000, EI=2*10^7, A=0.01, rho=7600:

> y := x -> piecewise(
        x<=15/2, x^3/120000 - 383/480000*x,
        x<=17/2, -x^4/96000 - 135/4096 + 16109/960000*x - 9/2560*x^2 + 77/240000*x^3,
        17/2<x, 257/12000 - 1057/120000*x + x^2/1000 - x^3/30000
);

y := proc (x) options operator, arrow; piecewise(x <= 15/2, 1/120000*x^3-383/480000*x, x <= 17/2, -1/96000*x^4-135/4096+16109/960000*x-9/2560*x^2+77/240000*x^3, 17/2 < x, 257/12000-1057/120000*x+1/100...y := proc (x) options operator, arrow; piecewise(x <= 15/2, 1/120000*x^3-383/480000*x, x <= 17/2, -1/96000*x^4-135/4096+16109/960000*x-9/2560*x^2+77/240000*x^3, 17/2 < x, 257/12000-1057/120000*x+1/100...

ローディングブリッジの振動を表す 偏微分方程式です。

> deq := rho*A*diff(w(x,t),t,t) + EI*diff(w(x,t),x,x,x,x)=0;

deq := rho*A*(diff(w(x, t), `$`(t, 2)))+EI*(diff(w(x, t), `$`(x, 4))) = 0

はじめに、pdsolve コマンドを使い、偏微分方程式を解いてみます。

> sol:= pdsolve(deq,w(x,t));

sol := `&where`(w(x, t) = _F1(x)*_F2(t), [{diff(_F1(x), `$`(x, 4)) = _c[1]*_F1(x), diff(_F2(t), `$`(t, 2)) = -EI*_c[1]*_F2(t)/(rho*A)}])

2つの常微分方程式が得られました。この2つの常微分方程式の解の積が、PDE の解を与えます。x の微分方程式を選択し、deq1 に代入します。同様に、時間 t の微分方程式を選択し、deq2 に代入します。

> deq1:= op(select(has,op([2,1],sol),x));

deq1 := diff(_F1(x), `$`(x, 4)) = _c[1]*_F1(x)

> deq2:= op(select(has,op([2,1],sol),t));

deq2 := diff(_F2(t), `$`(t, 2)) = -EI*_c[1]*_F2(t)/(rho*A)

1番目の微分方程式の解

はじめに、x の微分方程式を解きます。橋は、端点の曲げモーメントのない2つのベアリングに乗っている棒と考える事ができます。橋の境界条件は、以下の通りです。

> boundary := y(0) = 0, y(L) = 0, D(D(y))(0) = 0, D(D(y))(L) = 0:

これらの条件から、Sturm-Liouvill の固有値問題を得ます。

> sol1:= rhs(dsolve(deq1,_F1(x)));

sol1 := _C1*exp(-I*_c[1]^(1/4)*x)+_C2*exp(-_c[1]^(1/4)*x)+_C3*exp(I*_c[1]^(1/4)*x)+_C4*exp(_c[1]^(1/4)*x)

_c1^(1/4) lambda と置き換え、三角関数に変換すると、次式を得ます。

> sol1a:= convert( subs(_c[1]^(1/4)=lambda, sol1 ), trig);

sol1a := _C1*(cos(lambda*x)-I*sin(lambda*x))+_C2*(cosh(lambda*x)-sinh(lambda*x))+_C3*(cos(lambda*x)+I*sin(lambda*x))+_C4*(cosh(lambda*x)+sinh(lambda*x))

定数 a, b, c, d を導入し、三角関数を整理します。

> sol1b := a*cos(lambda*x) + b*sin(lambda*x) + c*cosh(lambda*x) + d*sinh(lambda*x);

sol1b := a*cos(lambda*x)+b*sin(lambda*x)+c*cosh(lambda*x)+d*sinh(lambda*x)

一般解に境界条件を与えます。線形方程式のシステムが得られます。

> eq := [eval(subs(x=0, sol1b)), eval(subs(x=0, diff(sol1b, x, x))), eval(subs(x=L, sol1b)), eval(subs(x=L, diff(sol1b, x, x)))];

eq := [a+c, -a*lambda^2+c*lambda^2, a*cos(10*lambda)+b*sin(10*lambda)+c*cosh(10*lambda)+d*sinh(10*lambda), -a*cos(10*lambda)*lambda^2-b*sin(10*lambda)*lambda^2+c*cosh(10*lambda)*lambda^2+d*sinh(10*lam...eq := [a+c, -a*lambda^2+c*lambda^2, a*cos(10*lambda)+b*sin(10*lambda)+c*cosh(10*lambda)+d*sinh(10*lambda), -a*cos(10*lambda)*lambda^2-b*sin(10*lambda)*lambda^2+c*cosh(10*lambda)*lambda^2+d*sinh(10*lam...

> M:= linalg[genmatrix](eq,[a,b,c,d]);

M := matrix([[1, 0, 1, 0], [-lambda^2, 0, lambda^2, 0], [cos(10*lambda), sin(10*lambda), cosh(10*lambda), sinh(10*lambda)], [-cos(10*lambda)*lambda^2, -sin(10*lambda)*lambda^2, cosh(10*lambda)*lambda^...

微分方程式の解が 0 の場合には、興味がありません。そこで、線形連立方程式が1つ以上の解をもつように、lambda を決定します。行列式が 0 になる lambda を見つけます。

> tmp:=linalg[det](M);

tmp := -4*sin(10*lambda)*lambda^4*sinh(10*lambda)

tmp=0 を solve コマンドで解くと、1つの解 lambda = 0 (主値)しか返しません。しかしこの値は、解が定数になってしまうので、意味がありません。そこで sin 関数の根を、次のように置きます。

> assume(n,integer):
eigenvalue := lambda = (n*Pi)/L;

eigenvalue := lambda = 1/10*n*Pi

求まった固有値を、線形連立方程式に代入し、それを解きます。

> MM:= subs({eigenvalue}, convert(M, Matrix));

MM := Matrix([[1, 0, 1, 0], [-1/100*n^2*Pi^2, 0, 1/100*n^2*Pi^2, 0], [cos(n*Pi), sin(n*Pi), cosh(n*Pi), sinh(n*Pi)], [-1/100*cos(n*Pi)*n^2*Pi^2, -1/100*sin(n*Pi)*n^2*Pi^2, 1/100*cosh(n*Pi)*n^2*Pi^2, 1...

> linalg[linsolve](MM,[0,0,0,0]);

vector([0, _t[1], 0, 0])

ベクトル [ a, b, c, d ] と比較すると、b の値は任意で、a, c, d の値はゼロになります。各々の 固有値 に対する解を得ます。

> sol1c:=subs(a=0,c=0,d=0,sol1b);

sol1c := b*sin(lambda*x)

微分方程式の一般解は、任意定数 bn での、固有値の解の和になります。

> Sum(subs(eigenvalue, b=b[n], sol1c), n=1..infinity);

Sum(b[n]*sin(1/10*n*Pi*x), n = 1 .. infinity)

> u:= subs(eigenvalue, b=b[n], sol1c);

u := b[n]*sin(1/10*n*Pi*x)

2番目の微分方程式の解

2番目の微分方程式を解きます。

> sol2:= dsolve(deq2,_F2(t));

sol2 := _F2(t) = _C1*sin(EI^(1/2)*_c[1]^(1/2)*t/(rho^(1/2)*A^(1/2)))+_C2*cos(EI^(1/2)*_c[1]^(1/2)*t/(rho^(1/2)*A^(1/2)))

_c[1] lambda^4 を代入し、lambda に固有値関数 eigenvalue を代入します。

> v:= rhs( subs(_c[1]=lambda^4, eigenvalue, sol2));

v := _C1*sin(1/10000*EI^(1/2)*10000^(1/2)*(n^4*Pi^4)^(1/2)*t/(rho^(1/2)*A^(1/2)))+_C2*cos(1/10000*EI^(1/2)*10000^(1/2)*(n^4*Pi^4)^(1/2)*t/(rho^(1/2)*A^(1/2)))

解の結合

私たちの本当の問題の解は、2つの常微分方程式の解の積で求まります。

> w := unapply(Sum(u*v, n=1..infinity), x, t);

w := proc (x, t) options operator, arrow; Sum(b[n]*sin(1/10*n*Pi*x)*(_C1*sin(1/10000*EI^(1/2)*10000^(1/2)*(n^4*Pi^4)^(1/2)*t/(rho^(1/2)*A^(1/2)))+_C2*cos(1/10000*EI^(1/2)*10000^(1/2)*(n^4*Pi^4)^(1/2)*...

定数の値は、初期条件により求まります。変位の初速度はゼロなので、変位を微分し、t=0 で評価します。

> eval(subs(t=0, diff(w(x,t),t)));

Sum(1/10000*b[n]*sin(1/10*n*Pi*x)*_C1*EI^(1/2)*10000^(1/2)*(n^4*Pi^4)^(1/2)/(rho^(1/2)*A^(1/2)), n = 1 .. infinity)

これは、定数 _C1 がゼロであることを示します。変位の初期値は

定数 _C2 b[n] は、変位の初期値により変わります。変位の初期値 y をフーリエ級数として変換します。

フーリエ級数を計算

> f:= unapply( piecewise(x<0, -y(-x), y(x)), x);

f := proc (x) options operator, arrow; piecewise(x < 0, -piecewise(-x <= 15/2, -1/120000*x^3+383/480000*x, -x <= 17/2, -1/96000*x^4-135/4096-16109/960000*x-9/2560*x^2-77/240000*x^3, 17/2 < -x, 257/120...f := proc (x) options operator, arrow; piecewise(x < 0, -piecewise(-x <= 15/2, -1/120000*x^3+383/480000*x, -x <= 17/2, -1/96000*x^4-135/4096-16109/960000*x-9/2560*x^2-77/240000*x^3, 17/2 < -x, 257/120...f := proc (x) options operator, arrow; piecewise(x < 0, -piecewise(-x <= 15/2, -1/120000*x^3+383/480000*x, -x <= 17/2, -1/96000*x^4-135/4096-16109/960000*x-9/2560*x^2-77/240000*x^3, 17/2 < -x, 257/120...

y のフーリエ級数を求めます。y は奇関数なので、その係数 a[n] はゼロになります。

> a[0] := (2*int(f(x)*cos((0*2*Pi*x)/(2*L)), x=-L..L))/(2*L);

a[0] := 0

> a[n] := simplify((2*int(f(x)*cos((n*2*Pi*x)/(2*L)), x=-L..L))/(2*L));

a[n] := 0

> b[n] := simplify((2*int(f(x)*sin((n*2*Pi*x)/(2*L)), x=-L..L))/(2*L));

b[n] := 5*(-cos(3/4*n*Pi)+cos(17/20*n*Pi))/(n^5*Pi^5)

よってフーリエ級数は

> Sum(b[n]*sin((n*2*Pi*x)/(2*L)), n=1..infinity);

Sum(5*(-cos(3/4*n*Pi)+cos(17/20*n*Pi))*sin(1/10*n*Pi*x)/(n^5*Pi^5), n = 1 .. infinity)

最初の 10 項の和を求め、級数展開と元の関数を比較します。

> plot([sum(b[n]*sin((n*2*Pi*x)/(2*L)), n=1..10), y(x)], x=-L..L, color=[red, blue], thickness=2);

[Plot]

私たちの問題の解を求めるため、係数の値を、微分方程式の解に代入します。

> W:= subs(_C1=0, _C2=1, w(x,t));

W := Sum(5*(-cos(3/4*n*Pi)+cos(17/20*n*Pi))*sin(1/10*n*Pi*x)*cos(1/10000*EI^(1/2)*10000^(1/2)*(n^4*Pi^4)^(1/2)*t/(rho^(1/2)*A^(1/2)))/(n^5*Pi^5), n = 1 .. infinity)

解のグラフ表示

橋の振動を視覚化するため、解に定数の値を代入し、級数の5項目以降を切り捨てます。

> W:= unapply(value(subs(val,infinity=5,W)),(x,t));

W := proc (x, t) options operator, arrow; -0.3004704472e-2*sin(.3141592654*x)*cos(50.63000231*t)+0.3001161369e-3*sin(.6283185308*x)*cos(202.5200093*t)-0.5806273019e-4*sin(.9424777962*x)*cos(455.670020...W := proc (x, t) options operator, arrow; -0.3004704472e-2*sin(.3141592654*x)*cos(50.63000231*t)+0.3001161369e-3*sin(.6283185308*x)*cos(202.5200093*t)-0.5806273019e-4*sin(.9424777962*x)*cos(455.670020...

> plot3d(W(x,t), t=0..0.2, x=0..L, axes=boxed, numpoints=2000,
title="Oscillation of the Bridge",
labels=["Time [s]","Distance [m]","Displacement [m]"]);

[Plot]

結論

Maple を使い、荷重が 突然 離れたときの、ローディングブリッジの振動を計算できることを 見てきました。また、指定された振幅を超えるかどうかを、見ることもできます。

参考文献

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