025_BeamAnalysis.mw

梁の解析

イントロダクション

このアプリケーションは弾性梁の変位とその傾き、曲げモーメント、せん断力を解きます。梁解析プログラムは、一般的で、広範囲の支持条件に適用できます。また、ここでは MacCauley あるいは特異関数の原理を使用しています。このプログラムは Maple のプログラミング言語のパワーを示しています。

[Inserted Image]

> restart;

いくつかの有効なマクロ

以下はそれぞれの特異関数に通例の名前をつけ、境界条件の定義を単純化しています。

> dist_force := (x) -> Heaviside(x): # distributed force
conc_force := (x) -> Dirac(x): # concentrated force
conc_moment := (x) -> Dirac(1,x): # concentrated moment
D2 := (x) -> D(D(x)): # short form for 2nd derivative
D1 := (x) -> D(x): # short form for 1st derivative

梁解析プログラム

以下のサブルーチンは、弾性梁の変位とその傾き、曲げモーメント、せん断力を解きます。 アルゴリズムは、McCauleyの加圧関数を使った荷重条件、支持の境界条件を引数に、 4次の支配境界値問題を形式化したものです。 BVPは重ね合わせの方法を使って解析的に解きます。

入力:

MacCauley関数の概念を使った荷重条件のセット

梁の支持を指定する境界条件のセット

出力:

目標の解を含むテーブル

> beam_analyze := proc(loadset,bcset)
 local N,i,yy,th,m,v,q,de,outab;
 # Form differential equation
 de := EI*diff(y(x),x$4) = sum(loadset[i],i=1..nops(loadset)):

 dsolve({de} union bcset,y(x)): # Solve boundary value problem
 yy := rhs(%): # Extract deflection
 th := diff(yy,x): # Extract slope
 m := diff(yy,x$2): # Extract moment
 v := diff(yy,x$3): # Extract shear

 print(`Boundary Value Problem to Solve is`);
 print(lhs(de)=op(loadset),bcset);
 outab[`deflection`]:= simplify(yy):
 outab[`slope`] := simplify(th):
 outab[`moment`] := simplify(m):
 outab[`shear`] := simplify(v):
 outab:
end proc:

梁解析プログラムを使って

問題の記述

梁は単純に支持され、梁全体に一様に分布した荷重がかかっています。 荷重の分布は x = a から始まっており , x = b には力が集中し, x = cにはモーメントが集中しています。 beam-analyzeサブルーチンを使用し、この問題を解きます。

荷重と支持

荷重と支持 (単純支持) の定義

> loads := {-w*dist_force(x-a), -v*conc_force(x-b), m*conc_moment(x-c)}:
supports := {D2(y)(0) = 0, D2(y)(L)=0,y(L)=0,y(0)=0}:

パラメータの値

外部荷重の位置と大きさ

梁の力の分布は x = a から始まり、 次のように伝播します。

力は x = b に集中
モーメントは x = c に集中

分布荷重: w N/m

集中荷重: v N

集中モーメント: m N.m

材質の剛性: EI N.m^2

> L:=10: a:=0: b:=2: c:=4: v:=500: w:=500: m:=10000: EI:=10^6:

問題を解きます。

> sol := beam_analyze(loads,supports):

`Boundary Value Problem to Solve is`

1000000*(diff(y(x), `$`(x, 4))) = (-500*Heaviside(x), -500*Dirac(x-2), 10000*Dirac(1, x-4)), {y(0) = 0, `@@`(D, 2)(y)(0) = 0, y(10) = 0, `@@`(D, 2)(y)(10) = 0}

> newsol:=convert(sol,table):

ポスト処理

変位ダイアグラム

> plot(newsol[deflection], x=0..L, title=`Beam Deflection`, axes=BOXED);

[Plot]

せん断力ダイアグラム

> plot(newsol[shear], x=0..L, title=`Shear Force Diagram`, axes=BOXED);

[Plot]

曲げモーメント・ダイアグラム

> plot(newsol[moment], x=0..L, title=`Bending Moment Diagram`, axes=BOXED);

[Plot]

3D 剛性解析

剛性の関係式とパラメータ範囲を定義

> L:=10: a:=0: b:=2: c:=4: v:=500: w:=500: m:=10000: unassign('EI'):
sol := beam_analyze(loads,supports):

`Boundary Value Problem to Solve is`

EI*(diff(y(x), `$`(x, 4))) = (-500*Heaviside(x), -500*Dirac(x-2), 10000*Dirac(1, x-4)), {y(0) = 0, `@@`(D, 2)(y)(0) = 0, y(10) = 0, `@@`(D, 2)(y)(10) = 0}

> stiff:=convert(sol,table):

> EIrange := 1e5 .. 2e6;

EIrange := 0.1e6 .. 0.2e7

変位対剛性

> plot3d(stiff[deflection],x=0..10,EI=EIrange,axes=BOXED);

[Plot]

曲げモーメント対剛性

> plot3d(stiff[moment],x=0..10,EI=EIrange,axes=BOXED);

[Plot]

解析問題の例題

荷重分布が変動する場合の解析

以下の条件で位置に線形に依存するように外部荷重分布を変化させます。

分布荷重: LHS(左端)=200 そして RHS(右端)=1000
集中荷重: 5500 ( x=2 m にて)

集中モーメント: なし

変位およびせん断力ダイアグラムをプロット

> L:=10: a:=0: b:=2: c:=0: v:=5500: w:=200+80*x: m:=0: EI:=10^6:
sol:=beam_analyze(loads,supports):

`Boundary Value Problem to Solve is`

1000000*(diff(y(x), `$`(x, 4))) = (0, -(200+80*x)*Heaviside(x), -5500*Dirac(x-2)), {y(0) = 0, `@@`(D, 2)(y)(0) = 0, y(10) = 0, `@@`(D, 2)(y)(10) = 0}

> newsol:=convert(sol,table):

> plot(newsol[deflection],x=0..L,title=`Beam Shear`,axes=BOXED);

[Plot]

> plot(newsol[shear],x=0..L,title=`Beam Deflection`,axes=BOXED);

[Plot]

3D 集中荷重解析

集中荷重を 0 から 10000 まで変化させたときの変位の変化を観測します。分布荷重を取り除きます。すなわち w=0.

> newsol:=convert(eval(op(subs(L=10,a=0,b=2,c=0,w=0,m=0,EI=10^6,op(sol)))),table):

> L:=10: a:=0: b:=2: c:=0: w:=0: m:=0: EI:=10^6: unassign('v'):
sol:=beam_analyze(loads,supports):

`Boundary Value Problem to Solve is`

1000000*(diff(y(x), `$`(x, 4))) = (0, -v*Dirac(x-2)), {y(0) = 0, `@@`(D, 2)(y)(0) = 0, y(10) = 0, `@@`(D, 2)(y)(10) = 0}

> newsol:=convert(sol,table):
plot3d(newsol[deflection],x=0..10,v=0..10000,axes=BOXED);

[Plot]

設計問題の例題

1% の変位の許容範囲で分布した力の最大値を求めます。

以下の荷重条件が与えられたとします。

荷重分布 : x = 0 から x = 10まで

> defl:=convert(eval(op(subs(L=10,a=0,b=0,c=0,v=0,m=0,EI=10^6,op(sol)))),table):

> L:=10: a:=0: b:=0: c:=0: v:=0: m:=0: EI:=10^6: unassign('w'):
sol:=beam_analyze(loads,supports):

`Boundary Value Problem to Solve is`

1000000*(diff(y(x), `$`(x, 4))) = (0, -w*Heaviside(x)), {y(0) = 0, `@@`(D, 2)(y)(0) = 0, y(10) = 0, `@@`(D, 2)(y)(10) = 0}

> defl:=convert(sol,table):
maxdefl:=defl[slope]=0;

maxdefl := -1/6000000*w*(Heaviside(x)*x^3-15*x^2+250) = 0

最大変位はスパンの中央 x = 5 で発生します。

> newdefl:=simplify(subs(x=5,defl[deflection])=-0.01*L);

newdefl := -1/7680*w = -.1000000000

> maxforce:=eval(solve(newdefl,w)) * N;

maxforce := 768.*N

与えられた高さと荷重に対する最大の梁幅を探します。

以下の荷重条件が与えられたとします。

集中荷重 スパン中央において500 N
最大変位 スパン中央において長さの 1% が発生します。

高さ 0.1 m

1%の変位をサポートするアルミニウムの梁に対する、正方形の断面の面積を見つけます。

> cross:=convert(eval(op(subs(L=10,a=0,b=5,c=0,v=500,m=0,w=0,EI=Youngs*(base*height^3)/12,op(sol)))),table):

> L:=10: a:=0: b:=5: c:=0: v:=500: m:=0: w:=0: EI:=Youngs*(base*height^3)/12:
sol:=beam_analyze(loads,supports):

`Boundary Value Problem to Solve is`

1/12*Youngs*base*height^3*(diff(y(x), `$`(x, 4))) = (0, -500*Dirac(x-5)), {y(0) = 0, `@@`(D, 2)(y)(0) = 0, y(10) = 0, `@@`(D, 2)(y)(10) = 0}

> cross:=convert(sol,table):
maxhgt:=solve(simplify(subs(x=5,height=.1,Youngs=69E9,cross[deflection])=-.01*L),base);

maxhgt := 0.1811594202e-1

梁の境界条件を両端固定に変更し、先程と同じ解法を行います。

新しい一般解を決定するために、梁解析マクロを再実行します。
新しい支持条件は以下のように示されています。

> supports := {D(y)(0) = 0, D(y)(L)=0,y(L)=0,y(0)=0};

supports := {y(0) = 0, y(10) = 0, D(y)(0) = 0, D(y)(10) = 0}

> L:=10: a:=0: b:=5: c:=0: v:=500: m:=0: w:=0: EI:=Youngs*(base*height^3)/12:

> clamped := beam_analyze(loads,supports):

`Boundary Value Problem to Solve is`

1/12*Youngs*base*height^3*(diff(y(x), `$`(x, 4))) = (0, -500*Dirac(x-5)), {y(0) = 0, y(10) = 0, D(y)(0) = 0, D(y)(10) = 0}

問題を解きます。

> cross:=convert(clamped,table):

> maxhgt:=solve(simplify(subs(x=5,height=.1,Youngs=69E9,cross[deflection])=-.01*L),base);

maxhgt := 0.4528985506e-2